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ABSTRACT 

The  application  of  the  finite  difference-time  domain  (FD-TD)  technique  to  prob¬ 
lems  in  electromagnetic  interference  and  radiation  is  examined,  and  enhancements  to 
the  method  in  order  to  increase  its  accuracy  and  efficiency  are  considered.  The  pri¬ 
mary  motivation  for  studying  electromagnetic  radiation  from  various  structures  which 
compose  a  computer  is  that  computers  are  now  operating  at  faster  speeds  (i.e.  lOO’s  of 
MHz)  and  consuming  more  power  resulting  in  significantly  increased  levels  of  power 
at  higher  frequencies.  Three  electromagnetic  interference  problems  associated  with 
computers  are  analyzed.  First,  electromagnetic  radiation  from  VLSI  heatsink  con¬ 
figurations  is  examined.  Next,  electromagnetic  emissions  from  modules-on-backplane 
configurations  are  analyzed.  The  third  problem  involves  the  analysis  of  coupling  of 
electromagnetic  energy  by  wires  exiting  from  metallic  enclosures.  In  addition  to  these 
interference  problems,  the  effects  of  non-uniform  substrates  on  radiation  by  microstrip 
patch  antennas  are  studied.  In  particular  the  use  of  non-uniform  substrates  to  reduce 
loss  due  to  surface  waves  is  examined. 

The  electromagnetic  radiation  from  a  VLSI  chip  package  and  heatsink  structure 
is  analyzed  using  the  finite  difference-time  domain  technique.  The  dimensions  of  a 
typical  configuration  calls  for  a  multi-zone  gridding  scheme  in  the  FD-TD  algorithm 
to  accommodate  fine  grid  cells  in  the  vicinity  of  the  heatsink  and  package  cavity  and 
sparse  gridding  in  the  remainder  of  the  computational  domain.  The  issues  pertaining 
to  the  effects  of  the  heatsink  in  influencing  the  overall  radiating  capacity  of  the  con¬ 
figuration  are  addressed.  Analyses  are  facilitated  by  using  simplified  heatsink  models 
and  by  using  dipole  elements  as  sources  of  electromagnetic  energy  to  model  the  VLSI 
chip.  The  potential  for  enhancement  of  spurious  emissions  by  the  heatsink  structure 
is  illustrated.  For  heatsinks  of  typical  dimensions,  resonance  is  possible  within  the 
low  gigahertz  frequency  range.  The  exploitation  of  the  heatsink  as  an  emissions  shield 
by  appropriate  implementation  schemes  is  discussed  and  evaluated. 

The  electromagnetic  radiation  from  a  modules-on-backplane  configuration  found 
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in  computers  is  analyzed  by  means  of  the  finite-difference  time-domain  technique.  The 
issues  pertaining  to  the  effects  of  the  modules  in  influencing  the  overall  radiating  ca¬ 
pacity  of  the  configuration  are  addressed.  Analyses  are  facilitated  by  using  simplified 
modules-on-backplane  models  and  by  using  a  voltage  source  as  the  source  of  electro¬ 
magnetic  energy  at  a  connector  between  a  module  and  the  backplane.  The  potential 
for  enhancement  of  spurious  emissions  by  the  module-on-backplane  configuration  is 
demonstrated.  For  modules  of  typical  dimensions,  resonance  is  observed  at  the  hun¬ 
dreds  of  megahertz  frequency  range.  Comparisons  between  numerical  predictions  and 
measurements  are  conducted  to  validate  the  numerical  method. 

The  finite  difference-time  domain  technique  is  used  to  analyze  coupling  of  elec¬ 
tromagnetic  energy  through  metallic  enclosures  by  wires  and  cables  which  exit  them. 
Simplified  configurations  are  considered  in  order  to  focus  on  the  coupling  of  elec¬ 
tromagnetic  energy  by  a  wire  which  penetrates  the  metallic  enclosure.  A  lumped 
element  model  is  introduced  to  reduce  the  computational  domain  size  by  replacing 
the  external  wire  with  a  lumped  element  at  the  aperture  in  the  metallic  enclosure. 
The  lumped  element  consists  of  a  parallel  set  of  series  RLC  circuits,  and  is  designed 
to  match  the  radiation  admittance  of  the  external  wire  as  closely  as  possible.  The 
potential  for  significant  coupling  of  electromagnetic  energy  is  demonstrated  in  the 
hundreds  of  megahertz  frequency  range.  Comparisons  between  FD-TD  simulations 
using  the  lumped  element  model  and  FD-TD  simulations  modeling  the  the  entire  ex¬ 
ternal  wire  region  are  conducted  to  validate  the  numerical  method  and  evaluate  the 
lumped  element  model. 

The  analysis  of  microstrip  patches  on  non-uniform  substrates  with  special  em¬ 
phasis  on  the  reduction  of  surface  wave  losses  is  conducted  using  the  finite  difference¬ 
time  domain  technique.  The  alteration  to  the  uniform  substrate  involves  removing 
a  rectangular  ring  of  the  dielectric  substrate  which  surrounds  the  radiating  element. 
A  perfectly  matched  layer  absorbing  boundary  condition  is  used  for  enhanced  perfor¬ 
mance  especially  where  the  dielectric  substrate  intersects  the  outer  boundary.  Sim¬ 
ple  dipole  antennas  and  a  microstrip  patch  antenna  on  non-uniform  substrates  are 
examined.  It  is  demonstrated  that  the  guided  powers  for  these  configurations  are 
significantly  less  than  those  for  the  equivalent  uniform  substrate  configurations. 

Thesis  Supervisors:  Professor  J.  A.  Kong 
Dr.  Robert  T.  Shin 
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Chapter  1 
Introduction 


This  thesis  examines  the  application  of  the  finite  difference-time  domain  (FD-TD) 
technique  to  problems  in  electromagnetic  interference  and  radiation,  and  the  en¬ 
hancement  of  the  method  in  order  to  increase  its  accuracy  and/or  efficiency.  In  par¬ 
ticular,  three  electromagnetic  interference  problems  associated  with  computers  are 
analyzed.  First,  electromagnetic  radiation  from  VLSI  heatsink  configurations  is  ex¬ 
amined.  Next,  electromagnetic  emissions  from  modules-on-backplane  configurations 
is  analyzed.  The  third  problem  involves  the  analysis  of  coupling  of  electromagnetic 
energy  by  wires  exiting  from  metallic  enclosures.  In  addition  to  these  interference 
problems,  the  effects  of  non-uniform  substrates  on  microstrip  patch  antennas  is  stud¬ 
ied  with  an  emphasis  on  the  reduction  of  surface  wave  generation. 

1.1  Introduction  to  Finite  Difference- Time 
Domain  Technique 

The  finite  difference-time  domain  (FD-TD)  technique  was  introduced  by  K.  S.  Yee  in 
1966  [1],  and  has  since  been  applied  to  many  problems  in  electromagnetics  [2]- [83]. 
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This  method  directly  solves  Maxwell’s  time-dependent  curl  equations  by  discretizing 
them  in  space  and  time.  The  discretization  in  space  is  such  that  the  shortest  wave¬ 
lengths  of  interest  are  resolved  by  10  or  more  cells,  and  the  time  step  is  determined 
by  stability  criteria.  Electric  and  magnetic  fields  are  interlaced  on  a  spatial  grid,  and 
calculated  in  a  leap-frog  manner  which  allows  time  marching. 

The  FD-TD  method  has  grown  in  popularity  and  usefulness  as  computers  have 
grown  in  size  and  speed.  The  FD-TD  method  has  filled  a  niche  in  the  electromagnetics 
community  by  offering  an  alternative  to  other  numerical  methods,  such  as  method 
of  moments  [84]-[91],  various  high  frequency  techniques  [92]-[97],  and  finite  elements 
[11],  for  analyzing  electromagnetic  phenomena.  The  FD-TD  technique  has  notable 
advantages  over  each  of  these  other  methods.  The  FD-TD  method,  finite  element 
method  and  method  of  moments  are  all  essentially  exact  methods,  in  that  Maxwell’s 
equations  are  solved  numerically  and  one  of  the  most  significant  sources  of  error 
introduced  is  due  to  the  finite  discretization  of  the  problem.  The  FD-TD  method 
is  most  similar  to  the  finite  element  method  and  is  in  fact  a  subset  of  the  finite 
element  method.  The  basic  condition  imposed  on  finite  elements  to  convert  to  finite 
difference  methods,  is  the  assumption  that  the  fields  are  assumed  constant  over  a 
single  cell,  instead  of  having  some  pre-defined  spatial  dependence  over  a  single  cell, 
such  as  linear  or  piece-wise  sinusoidal  dependencies.  With  finite  element  methods  the 
computational  domain  is  broken  into  a  grid  which  typically  conforms  to  the  object 
under  study,  and  the  fields  are  approximated  by  basis  functions  associated  with  each 
cell.  With  FD-TD,  the  computational  domain  is  also  broken  into  a  grid,  which  more 
often  than  not,  does  not  conform  to  the  object  under  study,  and  the  fields  are  assumed 
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to  be  constant  within  single  cells.  Since  the  grid  frequently  does  not  conform  to  the 
object  under  study,  the  discretization  of  the  object  is  only  approximate,  unlike  typical 
finite  element  implementations.  In  particular,  with  the  most  popular  rectangular  grid 
in  FD-TD,  curved  surfaces  are  approximated  by  staircase  approximations.  Both  of 
the  methods  are  based  on  the  differential  form  of  Maxwell’s  equations  and  inherently 
have  unbounded  computational  domains.  Hence,  it  is  necessary  to  limit  the  size 
of  the  computational  domain  by  applying  an  artificial  absorbing  boundary  condition 
which  simulates  unbounded  space  at  the  outer  boundary  of  the  computational  domain 
[67]-[83]. 

The  FD-TD  method  differs  from  the  method  of  moments,  in  that  the  former 
is  based  on  the  differential  form  of  Maxwell’s  equations  in  the  time  domain  and  the 
latter  is  based  on  the  integral  form  of  Maxwell’s  equations  in  the  frequency  domain. 
One  particular  advantage  of  the  method  of  moments  is  that  since  an  integral  equation 
is  used,  the  problem  is  bounded,  and  hence  the  only  errors  introduced  are  due  to  the 
discretization  of  the  currents,  as  no  artificial  truncation  of  the  computational  domain 
is  required.  Another  advantage  of  the  method  of  moments  is  that  once  the  impedance 
matrix  has  been  calculated  from  the  integral  equation  and  then  inverted,  the  scattered 
fields  for  different  excitation  configurations  can  be  calculated  relatively  quickly. 

Certain  advantages  of  the  FD-TD  method  are  that  since  it  is  a  time  domain 
method,  it  can  be  used  to  solve  multiple  frequencies  simultaneously,  and  it  is  more 
flexible  in  implementing  modifications  to  models.  However,  one  drawback  of  excit¬ 
ing  multiple  frequencies  is  that  if  the  structure  is  highly  resonant  and  the  resonant 
frequency  is  excited  then  achieving  steady  state  will  require  significantly  more  time 


24 


CHAPTER  1.  INTRODUCTION 


steps.  The  number  of  time  steps  required  to  obtain  steady  state  will  be  propor¬ 
tional  to  the  Q  of  the  resonant  structure.  However,  extrapolation  schemes  based  on 
Prony’s  method  [114]  and  other  methods  [115]  have  been  developed  to  obtain  the  late 
time  response  of  highly  resonant  structures,  thus  reducing  the  computational  bur¬ 
den  associated  with  analyzing  resonant  structures  with  the  FD-TD  technique.  For 
comparison,  the  method  of  moments  is  more  efficient  at  obtaining  monostatic  radar 
cross  sections  for  multiple  observation  angles  at  a  single  frequency,  while  FD-TD  is 
more  efficient  at  obtaining  monostatic  radar  cross  sections  for  multiple  frequencies  at 
a  single  observation  angle. 

The  high-frequency  techniques,  which  include  physical  optics,  and  physical,  ge¬ 
ometrical,  and  uniform  theories  of  diffraction,  are  popular  for  studying  electrically 
large  objects.  They  are  usually  the  only  feasible  method  for  studying  objects  which 
are  over  20  wavelengths  in  size,  because  the  previous  numerical  exact  methods  are 
all  memory  and  computationally  intensive,  and  current  computers,  in  general,  can¬ 
not  handle  problems  of  this  size.  However,  they  are  inapplicable  to  low  frequency  or 
electrically  small  problems  and  are  approximate  methods,  in  that  they  do  not  include 
various  higher  order  effects,  such  as  multiple  interactions.  So,  in  contrast  with  the 
previous  methods,  including  the  FD-TD  method,  these  high-frequency  techniques  are 
approximate  instead  of  exact  solutions. 

One  class  of  problems  which  have  been  investigated  using  the  FD-TD  technique 
is  prediction  of  electromagnetic  scattering  and  radar  cross  section  [2]-[18],[98]-[101]. 
This  class  of  problems  included  scattering  from  canonical  shapes,  surface  features, 
and  actual  aircraft.  In  addition,  calculation  of  radiation  characteristics  of  wire,  mi- 
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crostrip  patch,  and  other  antennas  has  been  performed  [19], [20].  Various  microwave 
components,  such  as  vias,  striplines,  microwave  filters,  and  microstrip  lines  have  been 
studied  [21]-[25].  The  characteristic  impedance  and  effective  wavenumber  of  guiding 
structures  have  been  calculated  using  the  FD-TD  technique  [26]-[29].  The  FD-TD 
method  has  been  used  in  the  analysis  of  electromagnetic  interference  problems  asso¬ 
ciated  with  aircraft  and  computers  [30]-[39],[102]-[105].  The  technique  has  also  been 
used  to  study  scattering  from  periodic  media  [40]- [45].  Since  the  conception  of  FD- 
TD  method  for  electromagnetics  in  1966,  it  has  been  extensively  modified  to  increase 
its  range  of  applicability,  its  accuracy,  and  its  efficiency.  The  first  implementation  of 
the  FD-TD  technique,  which  is  still  perhaps  the  most  popular,  utilized  a  regular  rect¬ 
angular  lattice  [2], [6].  One  difficulty  with  rectangular  grids  is  that  curved  or  slanted 
surfaces  are  approximated  by  staircases.  In  certain  problems  the  staircase  approxi¬ 
mation  can  cause  significant  errors  due  to  unwanted  scattering  from  the  staircases. 
In  an  effort  to  address  this  shortcoming  of  rectangular  grids,  other  grids  have  been 
investigated,  such  as  triangular  grids  [46], [47],  a  grid  in  spherical  coordinates  [48], 
and  a  generalized  non-orthogonal  grid  [49]-[51].  However,  each  of  these  alternative 
grids  have  not  been  used  as  extensively,  primarily  because  the  grid  generation  is  a 
difficult  task  and  typically  requires  extensive  user  input  and  requires  significantly 
more  memory  storage  since  a  description  of  each  non-uniform  node  in  the  grid  must 
be  kept.  An  additional  problem  which  arises  with  non-uniform  grid  dimensions  is 
that  the  stability  criteria  for  the  time  step  size  is  determined  by  the  smallest  cell  in 
the  grid.  If  there  are  small  cells  anywhere  within  the  grid,  the  time  step  must  be 
reduced  to  accommodate  those  small  cells;  this  will  in  turn  require  more  time  steps 
to  be  executed  to  achieve  steady  state. 
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In  addition,  modifications  have  been  made  to  the  rectangular  grid  to  allow  sub- 
gridding  for  increased  resolution  [52]- [55],  and  local  deformations  to  the  grid  to  better 
model  curved  surfaces  [8], [51].  Subgridding  is  particularly  useful  in  cases  where  there 
are  smaller  but  important  features  in  the  object  under  study  which  must  be  carefully 
modeled.  In  these  cases  a  smaller  grid  is  used  only  in  the  region  required.  Hence, 
subgridding  potentially  allows  higher  resolution  where  needed  without  excessively  in¬ 
creasing  the  memory  required  or  the  number  of  computations  per  time  step.  One 
drawback  with  subgridding  is  that  the  time  step  must  be  chosen  to  accommodate  the 
smallest  cell,  which  imposes  additional  computational  costs  similar  to  those  incurred 
with  non-uniform  grids.  The  subgridding  scheme  incurs  some  additional  complexity 
to  the  FD-TD  algorithm.  Local  deformations  to  the  rectangular  grid  reduces  or  elim¬ 
inates  the  problem  of  staircase  approximations  to  curved  surfaces.  The  generation  of 
the  deformed  grid  typically  requires  some  user  input,  however  it  is  far  easier  to  gen¬ 
erate  than  a  triangular  grid.  While  those  cells  which  include  deformed  edges  require 
additional  information,  the  majority  of  the  cells  in  the  computational  domain  will  re¬ 
main  standard  rectangular  cells  which  are  very  efficient  in  memory  and  computational 
usage. 

Regardless  of  which  grid  is  used,  an  absorbing  boundary  condition  (ABC)  is 
needed  at  the  outer  boundary  of  the  computational  domain,  in  order  to  truncate 
the  computational  domain  and  to  simulate  unbounded  space  [67]-[83].  One  of  the 
earlier  forms  of  ABC  was  to  actually  model  a  layer  of  lossy  material  at  the  outer 
boundary  which  hopefully  absorbed  most  of  the  outgoing  waves  and  incurred  min¬ 
imal  reflections.  This  solution  was  somewhat  undesirable  because  the  outer  layer 
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tended  to  be  fairly  large,  and  sometimes  imposed  a  severe  computational  penalty. 
So  ABCs,  which  are  approximations  to  the  one-way  wave  equation,  were  developed 
and  displayed  fairly  good  characteristics  and  were  efficient  [67]-[77].  Essentially,  these 
one-way  wave  equations  factor  the  wave  equation  into  inward  and  outward  traveling 
waves,  and  the  square  root  associated  with  these  one-way  wave  equations  is  approx¬ 
imated  using  various  polynomial  expansions.  These  ABCs  typically  perform  better 
at  angles  near  normal  incidence,  although  certain  implementations  can  be  tailored  to 
have  minimized  reflections  at  user  specified  angles.  However,  significant  reflections 
occur  at  angles  near  grazing  to  the  outer  boundary.  These  ABCs  have  been  applied 
to  the  various  types  of  grids.  Recently,  a  perfectly  matched  layer  (PML)  ABC  has 
been  developed  which  inserts  a  reflectionless  lossy  layer  of  finite  thickness  at  the  outer 
boundary  of  the  computational  domain  [78]- [83].  This  reflectionless  layer  is  achieved 
by  modifying  Maxwell’s  equations  (non-physical)  in  the  perfectly  matched  layer.  The 
performance  of  this  PML  ABC  has  been  exceptionally  good  with  reasonable  compu¬ 
tational  and  memory  costs.  In  particular,  the  angular  performance  of  the  PML  ABC 
appears  to  be  better  than  most  of  the  one-way  wave  equation  ABCs. 

The  earliest  implementations  of  the  FD-TD  technique  dealt  primarily  with  per¬ 
fectly  conducting  bodies,  but  has  since  been  extended  to  dielectric,  conductive  and 
magnetic  media  [2],  and  more  recently  to  dispersive  [9], [10], [56], [57],  gyrotropic  [58], 
and  anisotropic  media  [21].  Modifications  to  the  FD-TD  method  for  thin  dielec¬ 
tric  and  conducting  layers  have  been  performed  in  order  to  avoid  decreasing  the  cell 
size  to  the  thickness  of  the  dielectric  layer  [59]-[65].  In  addition,  modeling  of  small 
apertures  (sub-cell  dimensions)  in  perfectly  conducting  screens  has  been  recently  im- 
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plemented  in  the  FD-TD  method  through  the  use  of  modified  or  “smart”  FD-TD 
cells  [15].  In  each  of  these  prior  cases,  reducing  the  cell  size  locally  would  not  in  itself 
be  prohibitively  difficult,  however  the  associated  reduction  in  the  time  step  size  to 
accommodate  the  small  cells  would  be  computationally  cumbersome. 
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The  finite  difference-time  domain  technique  is  highly  accurate,  and  can  generate  virtu¬ 
ally  exact  solutions  for  scattering  and  radiation  problems  [1]-[12].  The  approximations 
associated  with  the  finite  difference-time  domain  technique  include  the  use  of  differ¬ 
ence  equations  in  place  of  differential  equations  and  the  use  of  an  absorbing  boundary 
condition  which  simulates  unbounded  space.  Theoretically,  the  computational  domain 
is  unbounded,  consisting  of  the  infinite  region  surrounding  the  scatterer.  However, 
an  absorbing  boundary  condition  can  be  utilized  at  an  artificial  outer  boundary  to 
make  the  computational  domain  finite  and  to  simulate  the  unbounded  space  beyond. 
Within  the  finite  computational  domain.  Maxwell’s  equations  are  discretized  on  a 
grid.  In  order  to  obtain  accurate  results,  the  spacing  between  the  nodes  in  the  grid 
must  be  on  the  order  of  a  tenth  of  a  wavelength.  Given  the  excitation  sources,  the 
electromagnetic  scattering  from  the  object  can  be  found  by  marching  in  time. 

The  most  popular  formulation  of  the  FD-TD  technique  is  Yee’s  algorithm  [1]. 
Yee  discretizes  the  electric  and  magnetic  fields  on  a  rectangular  grid  [See  Figure  1.1]. 
The  spatial  and  temporal  partial  derivatives  in  Maxwell’s  equations  are  replaced  with 
centered  difference  approximations.  Two  major  advantages  associated  with  Yee’s 
algorithm  are  that  objects  of  arbitrary  geometry  can  be  easily  modeled,  and  that 
the  algorithm  is  easily  implemented.  This  algorithm  has  second-order  accuracy  in 
both  time  and  space.  The  disadvantages  of  rectangular  grids  include  the  use  of 
staircase  approximations  for  curved  or  slanted  surfaces  and  the  use  of  a  rectangular 
outer  boundary.  Generally,  the  staircase  approximation  can  model  curved  surfaces 
reasonably  well.  However,  for  some  cases  such  as  when  strong  surface  waves  are 
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Figure  1.1:  Field  components  in  Yee’s  rectangular  grid. 

present,  a  very  fine  resolution  rectangular  grid  must  be  employed  to  obtain  accurate 
results. 

On  the  outer  boundary  of  this  grid,  an  absorbing  boundary  condition  (ABC) 
must  be  used  to  simulate  unbounded  space  [67]- [83].  Ideally,  the  absorbing  bound¬ 
ary  will  absorb  all  outgoing  waves,  thus  producing  no  reflections.  Realistically,  the 
absorbing  boundary  condition  is  chosen  to  minimize  the  reflections  over  all  possible 
angles.  The  absorbing  boundary  condition  must  be  well  posed  in  order  to  construct 
stable  difference  approximations.  Also  the  absorbing  boundary  condition  should  be 
local  in  time  and  space,  so  that  only  a  limited  number  of  time  steps  of  information 
near  the  outer  boundary  need  be  stored.  A  commonly  employed  local  absorbing 
boundary  condition  for  rectangular  grids  is  the  second-order  approximation  derived 
by  Engquist  and  Majda  [68].  In  treating  the  edges  of  the  outer  boundary,  a  first-order 
absorbing  boundary  condition  is  utilized  since  tangential  differences  are  unavailable. 
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This  combination  of  absorbing  boundary  conditions  has  been  shown  to  be  effective 
[2l,[67l,[77|. 

In  addition  to  the  boundary  conditions  at  the  outer  boundary,  boundary  con¬ 
ditions  at  material  interfaces  are  enforced.  In  particular,  at  perfectly  conducting 
surfaces  the  tangential  electric  fields  are  set  to  zero.  For  interfaces  between  different 
dielectric  and  magnetic  materials,  the  tangential  electric  and  magnetic  fields  are  kept 
continuous. 

In  using  the  finite  difference-time  domain  technique,  the  fields  propagate  through 
the  computational  domain  in  discrete  time  steps.  The  time  marching  is  executed  such 
that  the  electric  and  magnetic  fields  are  calculated  in  an  alternating  manner.  In  order 
to  insure  the  stability  of  the  time  stepping  algorithm,  the  size  of  the  time  step  must 
satisfy  Equation  1.1  for  three-dimensional  rectangular  grids  [6]. 

where  Ax,  Ay,  and  Az  are  the  dimensions  of  the  cell.  At  is  the  time  step,  and  c  is 
the  speed  of  light. 

Various  excitation  sources  such  as  plane  waves,  dipoles,  or  voltage  or  current 
sources  can  be  used.  The  time  dependence  of  these  sources  can  be  sinusoidal  or  a 
Gaussian  pulse.  The  Gaussian  pulse  can  be  used  for  analyzing  multiple  frequencies  by 
Fourier  transforming  the  tinie  response  of  the  system.  Frequently,  the  computational 
domain  will  be  separated  into  two  regions  when  a  plane  wave  excitation  is  used. 
In  the  inner  region  near  the  object,  the  total  fields  are  calculated,  whereas  in  the 
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outer  region,  the  scattered  fields  are  calculated.  At  the  boundary  between  the  two 
regions,  the  incident  field  is  added/subtracted  from  the  scattered/total  fields  needed 
to  calculate  the  total/scattered  fields  on  the  opposite  side  of  the  boundary. 

The  FD-TD  simulation  is  normally  begun  with  zero  fields  everywhere  within  the 
computational  domain.  At  each  following  time  step,  the  magnetic  fields  are  calculated 
in  terms  of  the  electric  and  magnetic  fields  of  the  previous  half  and  full  time  steps  using 
the  difference  equations  obtained  from  Maxwell’s  equations.  Next  the  electric  fields 
are  calculated  in  a  similar  manner.  The  absorbing  boundary  conditions  are  applied 
to  the  electric  fields  tangential  to  the  outer  boundary.  The  electric  fields  tangential 
to  any  perfectly  conducting  surfaces  are  set  to  zero.  The  excitation  source  can  be 
implemented  at  cells  containing  a  voltage  or  current  source,  or  a  dipole.  Alternatively 
with  a  plane  wave  source,  the  implementation  occurs  at  the  interface  between  the  total 
and  scattered  field  regions.  These  steps  are  essentially  repeated  until  steady  state  is 
reached  for  a  sinusoidal  excitation,  or  until  all  the  transient  scattered  fields  have 
propagated  out  of  the  computational  domain  for  a  Gaussian  pulse  excitation. 

During  the  FD-TD  simulation,  certain  electric  and  magnetic  fields  are  stored 
for  post-processing.  For  Gaussian  pulse  time-dependent  excitations,  it  is  necessary 
to  store  the  entire  time  response  so  that  it  can  be  Fourier  transformed  and  frequency 
domain  information  can  be  extracted.  With  sinusoidal  time-dependent  excitations, 
it  is  only  necessary  to  determine  the  magnitude  and  phase.  In  cases  where  the  near¬ 
field  response  at  a  particular  location  is  desired,  then  only  the  field  values  at  the  cell 
containing  that  location  need  be  stored.  Alternatively  with  voltage  or  current  sources, 
it  is  occasionally  desirable  to  store  the  voltage  and  current  at  the  source  location,  from 
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which  the  radiation  impedance  of  an  antenna  can  be  determined.  Frequently  the 
tangential  scattered  fields  on  a  surface  enclosing  the  object  of  interest  are  calculated. 
From  these  scattered  fields,  the  radiated  power  and  the  far-field  response  can  be  easily 
obtained  [106]-[109]. 
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1.3  Description  of  Thesis 

This  thesis  deals  exclusively  with  the  application  of  the  finite  difference-time  domain 
(FD-TD)  technique  to  solve  new  problems,  and  the  enhancement  of  the  method  in 
order  to  increase  its  accuracy  and/or  efficiency.  Following  the  introduction  in  Chap¬ 
ter  1,  the  finite  difference-time  domain  technique  is  described  in  detail  in  Chapter  2. 
Those  aspects  of  the  method  which  are  common  to  each  of  the  problems  examined 
are  addressed,  while  modifications  to  the  method  for  specific  problems  are  discussed 
in  the  chapters  describing  the  problem.  Difference  equations  for  free  space  and  di¬ 
electric/magnetic  media  are  derived,  as  well  as  difference  equations  for  the  treatment 
of  interfaces  between  dielectric/magnetic  media,  and  perfect  electric  and  magnetic 
conductors.  Stability  and  accuracy  issues  are  discussed.  In  addition,  the  implemen¬ 
tation  of  infinite  plane  wave  sources  as  well  as  local  voltage  or  current  sources  are 
examined.  The  treatment  of  the  outer  boundary  by  an  absorbing  boundary  condition 
is  also  reviewed.  The  conversion  of  near-field  data  to  far-field  data  is  also  addressed. 
Finally  some  limitations  of  the  method  and  its  overall  capabilities  are  discussed. 

The  next  three  chapters  concentrate  on  electromagnetic  interference  in  com¬ 
puter  systems.  The  extension  of  regulatory  electromagnetic  emissions  limits  to  the 
gigahertz  range  of  frequencies  for  high-performance  computing  equipment  necessitates 
caution  in  the  design  and  implementation  of  components  within  computing  systems 
[102], [104], [105].  The  avoidance  of  potentially  efficient  radiators  at  the  design  and 
development  phases  is  of  interest.  The  objective  is  to  minimize  the  need  for  cost- 
prohibitive  corrective  measures  commonly  invoked  when  problems  are  encountered  in 
tests  on  assembled  systems. 
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In  Chapter  3,  the  effects  of  a  heatsink  placed  over  a  chip  package  on  the  electro¬ 
magnetic  emissions  characteristics  of  the  configuration  are  examined  [31]- [33].  Such 
an  interest  is  justified  for  the  following  reasons.  The  heatsink  is  typically  metal¬ 
lic  and  in  close  proximity  to  the  chip  which  utilizes  large  amounts  of  power.  The 
heatsink  will  display  radiation  characteristics  similar  to  those  of  a  microstrip  patch 
antenna  [111]-[112].  The  electrical  dimensions  of  the  heatsink  structure  are  compa¬ 
rable  to  the  wavelengths  of  concern.  In  light  of  the  ever-increasing  power  levels  used 
in  high-performance  integrated  circuits,  the  heatsink  is  and  will  continue  to  be  an 
indispensable  component  of  the  VLSI  chip  packaging  configuration.  The  objectives 
are  to  identify  and  quantify  the  effects  of  the  heatsink  on  the  radiating  properties  of 
the  package  structure;  to  understand  and  exploit  the  radiation  mechanisms;  and  to 
evaluate  viable  heatsink  implementation  schemes  for  minimizing  the  overall  radiating 
capacity. 

The  electromagnetic  radiation  properties  of  a  heatsink  and  integrated  circuit 
(IC)  package  configuration  are  analyzed  by  means  of  suitable  models.  A  typical  heat¬ 
sink/package  configuration  of  interest  is  shown  in  Figure  1.2.  The  heatsink  may  be 
modeled  as  a  perfectly  conducting  rectangular  slab  positioned  over  a  finite-size  di¬ 
electric  medium  representing  the  chip  package.  The  chip  package  is  supported  in  turn 
by  a  dielectric  layer  of  infinite  extent  over  an  infinite  ground  plane  which  models  the 
substrate  or  printed  wiring  board  (PWB)  with  (at  least)  one  reference  layer.  The 
electromagnetic  sources  used  to  model  the  active  chip  include  electric  and  magnetic 
dipoles  of  vertical  and  horizontal  orientations  positioned  on  the  package/substrate 
dielectric  interface.  The  choice  of  such  elements  circumvents  the  problem  of  deriving 
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Figure  1.2:  Cut  away  view  of  a  VLSI  package/heatsink  configuration. 

a  rigorous  model  for  the  VLSI  circuitry  as  an  electromagnetic  source  and  eliminates 
factors  which  may  otherwise  clutter  the  understanding  and  interpretation  of  the  phys¬ 
ical  processes  associated  with  the  spurious  radiation  study.  The  resulting  model  of 
an  electromagnetically-coupled  slab  structure  is  shown  in  Figure  1.3.  The  use  of  a 
model  with  reduced  complexity  for  both  source  and  heatsink  is  to  facilitate  the  in¬ 
terpretation  of  results  based  on  underlying  physics.  The  extension  to  other  heatsink 
configurations  (e.g.,  circular  or  finned)  and  source  models  may  be  accomplished  with 
the  numerical  solution  employed.  With  sufficiently  fine  grids,  circular  heatsinks  may 
be  adequately  modeled.  The  model  in  Figure  1.3  is  readily  modified  to  model  vari¬ 
ous  heatsink  implementation  schemes  proposed  for  emissions  reduction  consideration. 
Field  strengths,  and  radiated  powers  are  computed  for  quantitative  analysis  and  eval¬ 
uation. 


Certain  technical  issues  in  the  application  of  the  FD-TD  technique  to  this  prob- 
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Figure  1.3:  Cross-sectional  view  of  simplified  package/heatsink  model. 

lem  will  be  addressed.  The  use  of  a  multizone  gridding  scheme  to  accommodate 
fine  grid  cells  in  the  vicinity  of  the  heatsink  and  package  cavity  and  sparse  grids  in 
the  remainder  of  the  computational  domain  [10], [31], [52]- [55].  This  modification  al¬ 
lows  increased  resolution  and  accuracy  where  needed,  while  not  sacrificing  too  much 
computational  power  in  regions  where  the  increased  resolution  is  unnecessary.  The 
implementation  of  the  electric  and  magnetic  dipoles  of  vertical  and  horizontal  orienta¬ 
tions  in  the  FD-TD  scheme  requires  normalization.  Since  the  dipoles  are  implemented 
by  exciting  the  fields  at  a  single  node  in  the  rectangular  grid,  the  dimensions  of  the 
dipole,  and  hence  their  dipole  moments,  are  not  precisely  defined.  The  error  associ¬ 
ated  with  the  dipole  size  ambiguity  is  small,  although  it  does  increase  with  increasing 
frequency.  So  the  necessary  normalization  can  be  accomplished  by  performing  an  ad- 
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ditional  FD-TD  simulation  of  a  simplified  configuration  and  relating  those  results  with 
analytical  solutions  [10], [31],  although  the  error  associated  with  these  small  current 
sources  is  small  [15]. 

The  radiation  enhancement  for  various  dipole  sources  are  examined  and  quanti¬ 
fied.  The  effects  of  the  height  of  the  heatsink  cavity  and  the  dielectric  loading  of  the 
cavity  are  examined.  The  changes  to  the  radiation  characteristics  due  to  the  addi¬ 
tion  of  fins  to  the  heatsink  are  analyzed.  Two  modifications  for  reducing  the  overall 
radiation  from  the  heatsink  configuration  are  considered:  the  use  of  metallic  posts 
for  grounding  the  heatsink  to  the  ground  plane  and  the  application  of  a  conducting 
gasket  around  the  lip  of  the  heatsink/ground  plane  cavity.  Finally  some  comments 
are  made  regarding  the  overall  radiation  characteristics  of  heatsinks,  the  effectiveness 
of  the  two  radiation  reduction  measures  examined,  and  the  application  of  the  FD-TD 
method  to  this  problem. 

In  Chapter  4,  the  electromagnetic  emissions  characteristics  of  a  modules-on- 
backplane  configuration  found  in  computers  are  examined  [34], [35].  The  modules-on- 
backplane  is  of  interest  because  it  is  a  common  packaging  hierarchy.  Also  a  typical 
modules-on-backplane  configuration  has  reference  plane  dimensions  which  are  appre¬ 
ciable  electrically  at  frequencies  in  the  hundreds  of  megahertz,  and  the  structure,  as 
a  whole,  may  function  as  an  antenna.  The  objective  is  to  identify  and  quantify  the 
radiation  mechanisms  in  a  modules-on-backplane  configuration  and  to  validate  the 
numerical  model  with  measurements. 

A  typical  configuration  for  the  modules-on-backplane  is  shown  in  Figure  1.4. 
The  elements  of  the  packaging  hierarchy  of  interest  in  this  study  are  the  backplane. 
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Figure  1.4;  A  typical  modules-on-backplane  configuration. 

the  connector  and  the  module.  The  modules  and  backplane  are  printed  wiring  boards 
with  at  least  one  reference  layer.  The  connector  contains  numerous  wires  for  carrying 
signals  and  for  providing  power  to  the  various  reference  layers.  Due  to  the  complex¬ 
ity  of  the  modules-on-backplane  configuration,  it  is  necessary  to  simplify  the  problem 
to  some  reasonable  level  of  detail.  The  modules  and  backplanes  are  populated  with 
numerous  chips  and  contain  many  signal  paths.  However,  the  single  largest  conduc¬ 
tor  or  signal  path  in  them  is  typically  one  of  the  ground  or  reference  planes.  These 
reference  planes,  because  of  their  larger  dimensions,  tend  to  enhance  radiation  more 
at  lower  frequencies  than  the  other  elements  associated  with  the  modules/backplane. 
For  this  reason,  the  modules  and  backplane  will  be  modeled  as  perfectly  conducting 
plates.  The  source  considered  is  the  transient,  inductive  voltage  developed  across  the 
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Figure  1.5:  Simplified  model  of  the  modules-on-backplane  configuration. 

connector  between  the  module  and  the  backplane  [113].  This  source  is  implemented 
by  applying  a  voltage  source  between  the  module  and  the  backplane.  The  result¬ 
ing  simplified  model  of  the  modules-on-backplane  configuration  is  shown  in  Figure 
1.5.  The  simplified  modules-on-backplane  and  source  configuration  avoids  many  of 
the  complex  interactions  inherent  in  a  real  modules-on-backplane  configuration  which 
would  otherwise  clutter  or  mask  the  fundamental  radiation  mechanisms.  The  exten¬ 
sion  to  other  modules-on-backplane  configurations  (e.g.,  multiple  modules,  enclosed) 
and  source  models  may  be  accomplished  with  the  numerical  method  employed. 

In  this  problem,  the  source  modeling  in  the  FD-TD  algorithm  for  this  configura¬ 
tion  is  of  some  concern.  The  transient  voltage  source  is  implemented  by  imposing  the 
electric  field  values  in  the  vertical  direction  at  the  location  of  the  connector  [24], [25]. 
This  implementation  models  a  matched  source  coaxial  cable  transmission  line  excita¬ 
tion. 


In  applying  the  FD-TD  technique  to  certain  of  the  more  resonant  modules-on- 
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Figure  1.6:  Measurement  setup  for  modules-on-backplane  configurations. 


backplane  configurations,  a  large  number  of  time  steps  (on  the  order  of  thousands)  are 
required  in  order  to  achieve  steady  state.  An  extrapolation  scheme  based  on  Prony’s 
method  is  used  in  conjunction  with  a  relatively  short  FD-TD  time  response  to  effi¬ 
ciently  obtain  the  late  time  response  [114].  This  extrapolation  method  was  applied 
to  the  time  response  of  the  current  and  voltage  on  the  connector.  The  FD-TD  simu¬ 
lation  is  executed  long  enough  such  that  the  remaining  time  response  approximately 
contains  only  a  few  decaying  complex  exponentials.  Prony’s  method  approximately 
determines  the  complex  amplitudes  and  complex  frequencies  of  these  exponentials. 
Once  the  complex  amplitudes  and  frequencies  are  determined,  then  the  time  response 
is  extrapolated  to  the  desired  extent. 

Measurements  were  performed  using  an  HP  8510  Network  Analyzer.  The  net¬ 
work  analyzer  was  connected  to  the  modules-on-backplane  configuration  to  be  tested 
by  a  test  port  cable  and  semi-rigid  coaxial  cable  [see  Figure  1.6].  In  each  of  the 
measurements,  the  center  conductor  of  the  coaxial  cable  was  soldered  to  the  module 
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under  test  such  that  the  coaxial  cable  was  in  contact  with  at  least  one  of  the  refer¬ 
ence  planes  in  the  module.  5-parameter  data  (5ii)  was  recorded  using  this  set-up 
and  are  used  to  derive  the  radiation  impedances  and  radiated  powers  of  the  various 
configurations  tested. 

Calculations  for  various  single  module-on-backplane  configurations  are  conducted 
and  compared  with  measurements  in  order  to  determine  their  correlation.  Next  mul¬ 
tiple  module-on-backplane  configurations  are  examined  and  compared  with  measure¬ 
ments.  A  partially  enclosed  configuration  is  analyzed,  and  FD-TD  predictions  and 
measurements  are  compared.  Finally  observations  are  made  about  the  overall  accu¬ 
racy  of  the  FD-TD  technique  and  the  general  radiation  characteristics  of  modules- 
on-backplane  configurations. 

In  Chapter  5,  the  coupling  of  electromagnetic  energy  by  wires  exiting  from 
metallic  enclosures  is  examined.  Metallic  enclosures  are  commonly  used  to  house 
computer  systems  and  usually  have  cables  and  wires  exiting  from  them.  A  simplified 
configuration  is  shown  in  Figure  1.7.  Since  the  metallic  enclosure  functions  both  as 
a  electromagnetic  shield  and  as  a  resonator  structure,  the  most  important  sources  of 
radiation  from  the  box  come  from  any  holes  in  the  box  used  for  ventilation  and  from 
any  cables  exiting  from  the  box  used  for  power  or  data  transfer  [15], [39].  In  particular 
for  cables,  the  currents  on  the  interior  surface  of  the  enclosure  can  couple  to  wires 
exiting  the  enclosure,  and  the  wire  can  function  as  an  antenna  in  the  hundreds  of 
MHz  frequency  range  [20], [66].  For  these  reasons,  the  coupling  of  electromagnetic 
energy  out  of  the  metallic  enclosure  by  cables  is  of  great  interest.  One  of  the  goals 
of  this  study  is  to  determine  the  efficiency  of  the  coupling  of  electromagnetic  energy 
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MODULES-ON-BACKPLANE  AS  SOURCE 


Figure  1.7:  A  simple  wire  exiting  a  metallic  enclosure  configuration. 

from  sources  inside  the  enclosure  to  the  region  outside  through  a  cable  exiting  the 
enclosure. 

In  analyzing  the  problem,  the  following  assumptions  and  simplifications  are 
used.  The  radiation  admittance  of  the  outer  portion  of  the  wire  is  assumed  to  be 
simple  enough  such  that  it  can  be  adequately  modeled  by  a  set  of  of  series  RLC 
circuits  which  are  in  parallel  [116]- [118].  Essentially,  one  RLC  circuit  is  used  to 
match  each  resonance  present  in  the  wire’s  radiation  admittance.  Hence,  an  important 
aspect  of  this  research  is  the  development  of  a  model  for  this  circuit  within  the  FD- 
TD  algorithm  [22], [23].  An  additional  assumption  is  that  the  radiation  mechanism 
is  primarily  due  to  the  current  on  the  wire  and  not  the  fields  traveling  through  the 
aperture,  since  the  aperture  size  is  only  being  modeled  in  a  simple  manner  and  will 

depend  on  the  FD-TD  cell  size.  The  enclosure  and  wires  are  all  assumed  to  be 
perfectly  conducting. 
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The  motivation  for  using  the  lumped  element  model  within  the  FD-TD  tech¬ 
nique  instead  of  a  full  FD-TD  simulation  is  to  obtain  significant  computational  sav¬ 
ings.  In  particular,  the  lumped  element  model  replaces  the  model  of  the  region  outside 
the  enclosure,  thus  eliminating  the  need  to  include  the  outer  region  in  the  FD-TD 
simulation.  The  computational  requirements  of  the  lumped  element  model  are  a 
mere  fraction  of  the  remaining  FD-TD  requirements,  on  the  order  of  tens  of  storage 
variables  compared  to  tens  of  thousands  of  storage  variables,  and  similar  ratios  for 
arithmetic  operations.  Hence,  if  the  computational  domain  can  be  reduced  by  50 
percent,  then  the  computational  savings  will  be  very  close  to  50  percent. 

In  addition,  if  the  external  region  is  simple  enough,  the  radiation  admittance  of 
the  wire  can  be  approximated  using  analytical  solutions  for  wires  over  infinite  ground 
planes.  Hence,  the  external  problem  can  be  solved  first,  and  then  the  internal  region 
can  be  analyzed  using  the  lumped  element  model  with  the  FD-TD  technique.  Also  by 
splitting  the  problem  into  two  regions,  the  probable  problem  of  a  significant  difference 
in  scale  between  features  in  the  two  regions  is  naturally  resolved,  as  each  region  can 
be  examined  separately  using  whatever  scale  is  most  appropriate  for  each  region. 

Two  simplified  configurations  are  considered  in  analyzing  the  real  configuration. 
The  simplest  model  involves  a  wire  over  an  infinite  ground  plane  being  excited  by 
a  coaxial  line.  This  first  model  is  used  to  validate  the  lumped  element  model,  since 
analytical  solutions  exist  for  a  complex  load  terminating  a  coaxial  transmission  line. 
This  model  will  also  provide  a  measure  for  how  well  the  lumped  element  can  model 
the  true  radiation  admittance.  The  next  model  involves  a  wire  connected  to  an 
infinite  metal  screen  exiting  through  a  small  aperture,  with  a  dipole  source  on  the 


1.3.  DESCRIPTION  OF  THESIS 


45 


“inner”  side  of  the  metal  screen.  In  each  of  these  two  cases,  the  radiation  admittance 
and  total  radiated  power  are  examined  for  FD-TD  simulations  modeling  the  entire 
structure  and  FD-TD  simulations  using  the  lumped  element  model.  Finally,  the 
usefulness  of  the  lumped  element  model,  and  the  overall  radiation  characteristics  of 
this  configuration  are  summarized. 

In  Chapter  6,  the  effects  of  a  non-uniform  substrate  on  the  radiation  character¬ 
istics  of  microstrip  patch  antennas  is  examined  using  the  finite  difference-time  domain 
technique  [147]-[152],[111],[119]-[144].  The  substrate  is  altered  by  the  removal  of  a 
rectangular  ring  from  the  dielectric  substrate  which  surrounds  the  patch  antenna  as 
shown  in  Figure  1.8.  A  microstrip  patch  antenna  is  examined  because  it  is  a  popular 
antenna  which  is  printed  on  a  substrate.  In  designing  microstrip  antennas,  some  fea¬ 
tures,  which  are  commonly  optimized,  include  minimizing  the  overall  antenna  size, 
maximizing  the  antenna’s  bandwidth  in  frequency  and  minimizing  the  losses  due  to 
surface  waves.  A  higher  dielectric  constant  substrate  can  be  used  to  reduce  the  an¬ 
tenna  size  as  the  wavelength  is  reduced  in  the  dielectric  with  increasing  dielectric 
constant  [111], [140], [141].  The  bandwidth  is  proportional  to  the  height  of  the  mi¬ 
crostrip  patch.  However,  in  achieving  these  two  goals,  the  loss  due  to  surface  waves 
is  increased  which  is  undesirable.  The  introduction  of  this  discontinuity  into  the  sub¬ 
strate  allows  a  significant  portion  of  the  power  in  the  surface  waves  to  be  reflected  at 
the  free  space/dielectric  interface  and  be  radiated  into  the  free  space  region. 

The  generation  of  surface  waves  is  undesirable  for  the  following  reasons.  First, 
power  coupled  into  guided  waves  is  power  which  is  not  radiated  which  lowers  the 
overall  radiation  efficiency  of  the  antenna  [140].  Also,  with  a  finite  size  substrate,  the 
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Figure  1.8:  Microstrip  patch  antenna  on  a  non-uniform  dielectric  substrate. 
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guided  waves  could  cause  unwanted  high  side  lobes  or  back  lobes  due  to  interaction 
with  the  truncation  of  the  substrate.  When  phased  arrays  are  mounted  on  a  common 
substrate,  significant  coupling  between  elements  can  occur  which  can  cause  limitations 
in  the  scanning  range  of  the  array  [147], [148]. 

A  perfectly  matched  layer  (PML)  absorbing  boundary  condition  (ABC)  [78]- 
[83]  is  utilized  instead  of  the  standard  one-way  wave  equation  absorbing  boundary 
condition  [67]-[74].  The  PML  ABC  is  chosen  for  its  superior  performance  in  general 
and  in  particular  for  its  handling  of  dielectric  interfaces  at  the  outer  boundary. 

In  comparing  microstrip  patch  antennas  over  uniform  and  non-uniform  sub¬ 
strates,  the  following  radiation  characteristics  were  examined.  The  effects  of  the 
addition  of  the  non-uniformity  on  resonant  frequency  and  bandwidth  are  observed. 
The  guided  power  for  each  configuration  is  also  calculated.  The  guided  power  is  ap¬ 
proximated  by  integrating  the  Poynting  power  passing  through  the  portion  of  a  closed 
surface  enclosing  the  antenna  which  also  lies  within  the  substrate.  The  non-uniformity 
is  shown  to  reduce  the  amount  of  power  coupled  into  surface  waves.  In  addition,  the 
gain  patterns  of  the  two  configurations  are  examined.  Finally,  the  advantages  and 
disadvantages  of  using  a  non-uniform  substrate  is  discussed. 

Chapter  7  summarizes  the  work  presented,  highlights  some  of  the  more  signifi¬ 
cant  results,  and  suggests  future  directions  for  research. 


CHAPTER  1.  INTRODUCTION 


Chapter  2 


Finite  Difference-Time  Domain 
Technique 

2.1  Introduction 


The  first  step  in  applying  the  FD-TD  technique  involves  approximating  Maxwell’s 
equations  in  differential  form  by  center  differences  in  space  and  time.  The  locations 
at  which  electric  and  magnetic  fields  are  calculated  are  positioned  on  some  sort  of 
grid.  For  three-dimensional  problems,  the  most  popular  grid  is  a  rectangular  mesh  [1]. 
A  significant  advantage  of  using  rectangular  grids  over  other  grids  [48]-[51]  is  greater 
simplicity.  However,  since  objects  are  discretized  on  a  rectangular  grid,  curves  and 
slanted  lines  are  approximated  by  staircases. 

To  begin  the  FD-TD  simulation,  all  the  fields  within  the  computational  domain 
are  initially  set  to  zero.  At  each  time  step,  the  electric  fields  are  calculated  in  terms 
of  the  electric  and  magnetic  fields  of  the  previous  time  step  using  the  difference 
equations  obtained  from  Maxwell’s  equations.  Next,  the  magnetic  fields  are  calculated 
in  a  similar  manner.  Boundary  conditions  are  enforced  at  the  outer  boundary  of 
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the  computational  domain  and  at  all  dielectric  and  conducting  interfaces.  At  the 
outer  boundary  a  second-order  absorbing  boundary  condition  is  utilized  in  order  to 
simulate  unbounded  space  beyond  the  computational  domain  [67]-[74].  The  tangential 
electric  fields  are  set  to  zero  at  the  conducting  surfaces,  and  the  tangential  electric 
and  magnetic  fields  are  kept  continuous  at  the  dielectric/ magnetic  boundaries.  The 
implementation  of  the  excitation  source  typically  requires  that  certain  electric  and 
magnetic  fields  be  updated  at  each  time  step.  These  steps  are  essentially  repeated 
until  steady  state  is  reached  for  a  sinusoidal  excitation,  or  until  all  the  transient 
scattered  fields  have  propagated  out  of  the  computational  domain  for  a  Gaussian 
pulse  excitation. 


2.2  Maxwell’s  Equations  in  Rectangular  Coordi¬ 
nates 


In  implementing  the  finite  difference-time  domain  technique,  Maxwell’s  equations 
must  be  discretized  in  space  and  time.  In  this  case,  Yee’s  lattice,  which  is  a  rectangular 
grid,  is  utilized  because  of  its  simplicity  [1].  Maxwell’s  equations  in  vector  differential 
form  in  an  isotropic  and  homogeneous  dielectric  and  magnetic  material  are 


V  xE  = 

rrH 

dt 

(2.1) 

VxH  = 

dD 

(2.2) 

V-D  = 

p 

(2.3) 

V-B  = 

Pm 

(2.4) 
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and  the  constitutive  relations  are 


where  eo  is  the  free-space  permittivity,  8.854  x  10“^^  F/m  and  /xq  is  the  free-space 
permeability,  47r  x  10“^  H/m.  In  addition,  and  Hr  are  respectively  the  relative 
permittivity  and  permeability,  p  and  pm  are  electric  and  magnetic  charge,  and  a  and 
am  are  electric  and  magnetic  conductivities,  respectively.  The  magnetic  charge  and 
conductivity  are  simply  mathematical  constructs  in  order  to  facilitate  the  analysis 
of  lossy  magnetic  materials.  Maxwell’s  divergence  equations,  (2.3)  and  (2.4),  are 
satisfied  in  the  finite  difference  scheme  by  applying  appropriate  initial  and  boundary 
conditions.  The  actual  difference  equations  used  in  the  FD-TD  technique  are  based 
upon  Maxwell’s  curl  equations,  (2.1)  and  (2.2),  and  the  constitutive  relations,  (2.5) 
and  (2.6). 

In  a  rectangular  coordinate  system.  Equations  (2.1),  (2.2),  (2.5)  and  (2.6),  can 
be  written  as  the  following  set  of  scalar  equations. 
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(2.8) 

(2.9) 
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CrCo 


dt 


dH^ 

dz 


dH, 

dx 


-  (^Ey 


dt 


dj^ 

dx 


(2.11) 

(2.12) 


In  free  space,  where  €r  =  /^r  =  1  a  =  cr,„  =  0,  Equations  (2.7)  to  (2.12) 
simplify  to 
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(2.18) 

2.3  Center  Differences  and  the  Unit  Cell 


The  following  notation  for  any  function  of  time  and  space  will  be  used  in  the  finite 
difference  equations, 

f{iAx,jAy,kAz,nAt)  =  P{i,j,k)  (2.19) 


The  partial  derivatives  in  space  and  time,  within  Maxwell’s  equations,  are  ap- 
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Figure  2.1:  Field  locations  on  three-dimensional  rectangular  grid, 
proximated  using  center  differences, 

am  .  ./K  +  a;/2)-/k-a{/2) 

di  ~  A{  ^  ’ 

The  electric  and  magnetic  field  components  are  interlaced  in  time,  and  are  calculated 
in  a  leap-frog  manner  (i.e.,  first  the  electric  fields  are  calculated,  then  the  magnetic 
fields  are  calculated,  and  the  sequence  is  repeated).  The  electric  and  magnetic  field 
components  are  interlaced  spatially  a  half  grid  cell  apart.  The  locations  of  the  field 
components  on  a  unit  cell  of  the  rectangular  grid  are  shown  in  Figure  2.1. 
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2.4  Difference  Equations  in  Free  Space 


The  finite  difference  equations  in  free  space,  are  obtained  by  applying  center  differ¬ 
ences  to  Equations  (2.13)  to  (2.18). 
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Rearranging  terms, 
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where  770  is  the  free  space  impedance  and  Cq  is  the  speed  of  light  in  free  space. 


2.5  Difference  Equations  for  Dielectric  and  Mag¬ 
netic  Media 


The  finite  difference  equations  in  an  isotropic,  homogeneous,  dielectric  and  magnetic 
material  are  obtained  by  discretizing  Equations  (2.7)  to  (2.12).  The  derivatives  are 
replaced  by  center  differences  and  the  electric  and  magnetic  field  terms  involving 
electric  and  magnetic  conductivities  are  approximated  by  using  an  average  of  the 
field  values  at  a  half  time  step  before  and  after  the  desired  time. 

The  finite  difference  equations  approximations  for  Maxwell’s  equations  in  di¬ 
electric  and  magnetic  material  are 
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Rearranging  terms, 
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2.6  Treatment  of  Interfaces  Between  Two  Media 


The  interface  between  any  two  media  is  chosen  to  occur  at  integer  nodes  {i,j,k), 
hence  the  electric  fields  are  tangential  and  magnetic  fields  are  normal  to  the  surface 
of  the  interface.  At  an  interface  in  the  x-y  plane  between  two  arbitrary  media,  the 
finite  difference  equations  for  E^,  Ey  and  H,  must  be  treated  carefully  [See  Figures 
2.2  and  2.3].  Since  the  normal  component  of  the  magnetic  field  is  not  continuous 
across  the  interface,  Hz  must  be  calculated  at  points  just  above  and  just  below  the 


Figure  2.2:  Interface  on  the  x-y  plane  between  two  arbitrary  media  viewed  in  the  x-z 
plane. 

interface. 
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Figure  2.3:  Interface  on  the  x-y  plane  between  two  arbitrary  media  viewed  in  the  y-z 
plane. 

Using  the  integral  form  of  Equation  (2.2)  to  calculate  Ey  yields  [See  Figure  2.2] 


<f  H-dl 
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=  IL  (“f + + Us, 


Converting  Equation  (2.49)  to  a  finite  difference  equation  yields 
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Applying  the  integral  form  of  Equation  (2.2)  to  calculate  E^  results  in  [See 
Figure  2.3] 


Transforming  Equation  (2.51)  to  a  finite  difference  equation  yields 
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2.7  Treatment  of  Perfect  Electric  Conductor 

The  boundary  condition  at  perfect  electric  conductors  states  that  the  tangential  elec¬ 
tric  field  must  be  zero, 


nxE  =  0  (2.53) 

Since  media  interfaces  are  chosen  to  occur  at  integer  node  (i,  j,  k)  locations,  it  is  clear 
from  Figure  2.1  that  the  electric  fields  are  always  tangential  to  the  surface.  Thus, 
in  the  finite  difference-time  domain  scheme,  the  boundary  condition  at  the  perfect 
electric  conductor  can  be  satisfied  by  simply  setting  these  electric  fields  equal  to  zero 
at  each  time  step. 


2.8  Treatment  of  Perfect  Magnetic  Conductor 

At  a  perfect  magnetic  conductor  the  tangential  magnetic  field  must  be  zero, 

hxH^O  (2.54) 

This  boundary  condition  is  satisfied  when  the  tangential  electric  field  at  the  surface 
of  the  perfect  magnetic  conductor  is  calculated.  One  way  of  obtaining  a  difference 
equation  for  the  tangential  electric  field  (i.e.,  Ex)  is  to  simplify  the  original  problem 
[Figure  2.4]  by  using  image  theory  [Figure  2.5].  The  difference  equation  for  Ex  at 
the  interface  between  a  perfect  magnetic  conductor  and  a  dielectric  and  magnetic 
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Figure  2.4:  Interface  on  the  x-y  plane  between  a  media  and  a  perfect  magnetic  con¬ 
ductor  viewed  in  the  y-z  plane. 


material  is  given  by 


(2.55) 
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Figure  2.5:  Image  problem  of  interface  on  the  x-y  plane  between  a  media  and  a  perfect 
magnetic  conductor  viewed  in  the  y-z  plane. 
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2.9  Stability  and  Accuracy 


The  choices  of  Arc,  Ay,  Az,  and  At,  are  motivated  by  reasons  of  accuracy  and 
stability  [6], [69].  In  general,  to  obtain  accurate  results.  Ax,  Ay  and  Az  must  be 
a  small  fraction  (~  1/10)  of  the  smallest  wavelength  in  any  media  expected  in  the 
model  or  of  the  smallest  dimension  of  the  scatterer.  To  ensure  the  stability  of  this 
time  stepping  algorithm.  At  must  satisfy  the  following  condition  [6]. 

=  +  i  +  ’  (2.56) 

2.10  Implementation  of  Excitation  Source 


The  two  basic  types  of  sources,  which  are  plane  waves  and  current  sources,  are  treated 
diflFerently.  When  the  excitation  source  is  assumed  to  be  either  a  sinusoidal  or  Gaus¬ 
sian  pulse  plane  wave,  the  computational  domain  is  separated  into  two  regions  in 
order  to  facilitate  the  treatment  of  the  excitation  source  [Figures  2. 6-2. 7]. 

For  a  finite  scatterer,  within  the  inner  region,  total  fields  are  calculated,  while 
in  the  outer  region,  only  scattered  fields  are  calculated  [Figure  2.6].  The  scattered 

fields  are  defined  to  be  the  difference  between  the  total  fields  and  the  incident  plane 
wave. 


Hs{x,z,t)  =  Htix,z,t)-Hi{x,z,t)  (2.57) 

Es{x,z,t)  =  Et{x,z,t)  -Ei{x,z,t) 


(2.58) 
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Figure  2.6:  Absorbing  boundary,  computational  domain  and  closed  surface  on  which 
fields  are  sampled  with  a  finite  scatterer. 


z 


Figure  2.7:  Absorbing  boundary,  computational  domain  and  surface  on  which  fields 
are  sampled  with  an  infinite  ground  plane  geometry. 
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The  incident  fields  are  subtracted  from  field  quantities  just  within  the  inner 
regions,  when  they  are  used  to  calculate  field  quantities  just  beyond  the  inner  region. 
The  incident  fields  are  added  to  field  quantities  just  outside  the  inner  regions,  when 
they  are  used  to  calculate  field  quantities  just  within  the  inner  region  [67]. 

The  finite  difference  equations  for  the  fields  in  the  vicinity  of  an  x-y  plane 
interface  between  the  total  and  scattered  field  regions  are 

£,V'(i  +  Ij, k)  =  ££(1  +  li, k)+no^  Wfm  +  j  +  |, *:) 

-  i/r/5  (i  +  1,  _  I,  *)J  _  I  (i  +  1, 4  +  1)  (2.59) 

+  +  1,  fc  +  i)]  -  +  1,2,  -  1)} 

J  +  2’  J  +  2>1^)  +  +  |>  +  2) 

+  £r‘(i,2  +  1,*  +  1)]  -  H^'(i,i  +  1,  -  1)}  (2.60) 

+  hi  +  f.*=)  - -  hi  +  i <=)] 

B^i-'(ij,k  +  h  =  e;AU,k  +  l) 

+10^  +  hi.  *  +  5)  -  -  J.i.*  +  5)]  (2.61) 

+  5-  +  5)  -  »«'(i'i  -  5.  *:  +  5)] 

+  1,  *  +  1)  =  +  1,  f:  +  1) 

{b,”  (i.2  +  1,*:  +  1)  -  [Ky(i,j  +  h  k) 


(2.62) 
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-  +  1,  (c)] }  -  i|l  +  1,  a:  +  1)  -  £”  (i,  j,  k  +  1)] 

ftV"  (i + j,  j,  A + j)  =  ffjr"  (j + i,  j,  A + i) 

[E:.(i  +  hi.  t  +  1)  -  E;,(i,j,k+  1)]  (2.63) 

+  hj.k  +  1)  -  [K(i  +  li.k)  -  K,ii  +  k.i.  *:)]} 

+  hi  +  5, «:)  =  +  hi+  h  <=) 

+^^[EZ{i+hi  +  hk)-K(i+hi.l=)\  (2-64) 

+  i  -  Kihi  +  h  ^)] 

For  an  infinite  ground  plane  configuration,  the  computational  domain  is  basi¬ 
cally  divided  into  the  region  above  the  ground  plane  and  the  region  below  the  ground 
plane  [Figure  2.7].  For  plane  wave  illumination  from  the  upper  half  space,  total  fields 
are  calculated  in  the  lower  region,  while  in  the  upper  region,  only  scattered  fields  are 
calculated.  In  this  case,  the  scattered  fields  are  defined  to  be  the  total  fields  minus 
the  incident  plane  wave  and  the  plane  wave  that  would  have  been  reflected  from  a 
uniform  infinite  ground  plane, 

Hs{x,z,t)  =  Ht{x,z,t)  -  Hi{x,z,t)  -  Hr{x,z,t)  (2.65) 

Es{x,z,t)  =  Et{x,z,t)  -  Ei{x,z,t)  -  Er{x,z,t)  (2.66) 


The  removal  of  the  reflected  plane  wave  from  the  upper  half  space  is  essential  because 
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the  reflected  plane  wave  is  inflnite  in  size  and  could  not  otherwise  be  adequately  mod¬ 
eled  within  the  flnite  computational  domain.  For  plane  wave  illumination  from  the 
lower  half  space  the  total  and  scattered  fields  are  simply  calculated  in  the  opposite 
regions.  The  implementation  of  a  plane  wave  is  basically  identical  to  the  implementa¬ 
tion  for  a  finite  scatterer  except  with  the  slightly  dilferent  definition  of  the  scattered 
field. 

The  implementation  of  current  sources  is  straightforward.  In  this  situation, 
total  fields  are  calculated  within  the  entire  computational  domain.  Both  electric  and 
magnetic  current  sources  can  be  implemented.  The  current  sources  are  implemented 
by  adding  the  appropriate  current  term  in  Maxwell’s  curl  equations. 
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where  J  and  Jm  are  the  electric  and  magnetic  current  densities.  Applying  center 


72 


CHAPTER  2.  FINITE  DIFFERENCE-TIME  DOMAIN  TECHNIQUE 


differences  results  in 

£?+'(*  +  i,},k)=  +lj,k)-  Ar%  J, 

+>?o|;  [nThi  +\,i+'rk)-  uTHi  +  5,  J  -  5.  *^)] 
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+h3  +  h  =  Hr  Hi +hj+h^)-  — 

^0 

+  i.3  + 1.*)  -  E2{i  +  li. k)] 


(2.78) 


For  a  current  source  in  the  z  direction,  the  electric  or  magnetic  current  density,  or 
Jmz,  is  assumed  to  have  length  Az  and  a  cross  section  of  Ax  Ay  [15]. 
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2.11  Absorbing  Boundary  Conditions 


An  absorbing  boundary  condition  (ABC)  at  the  outer  boundary  is  needed  to  make 
the  computational  domain  finite,  and  to  simulate  unbounded  space  beyond  the  com¬ 
putational  domain  [Figures  2.6-2. 7]  [67]-[74].  The  absorbing  boundary  condition  used 
for  this  problem  is  the  second-order  approximation  derived  by  Engquist  and  Majda 
[68], 


■  ^2  1/52  \ 

dr^  2  \dT?'^  dT^ ) 


w  =  0 


(2.79) 


where  it;  is  a  field  quantity  which  is  tangential  to  the  absorbing  boundary,  h  is  the 
normal  direction,  Ti,  T2  are  the  tangential  directions,  and  r  is  time  normalized  with 
respect  to  the  speed  of  light.  This  second-order  absorbing  boundary  condition  is 
applied  at  the  faces  of  the  computational  domain  [Figure  2.8].  The  second-order 
absorbing  boundary  condition  works  very  well  for  waves  which  are  at  or  near  normal 
incidence,  and  does  not  work  as  well  for  waves  which  are  incident  at  grazing  angles. 

In  treating  points  at  the  edges  of  the  outer  boundary,  a  first-order  absorbing 
boundary  condition  is  used  for  field  components  which  are  tangential  to  the  faces 
on  which  the  edges  reside  [Figure  2.8].  The  first-order  ABC  is  characterized  by  the 
following  equation, 


w  =  0 


(2.80) 


where  n  is  the  normal  direction. 
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Figure  2.8:  First-  and  second-order  absorbing  boundary  conditions  for  a  three- 
0  dimensional  computational  domain. 


^  The  absorbing  boundary  condition  at  the  edges  is  not  as  important  as  the  one  at 

the  faces,  because  there  are  fewer  points  on  the  edges  than  the  faces.  Thus,  in  general, 
using  the  first-order  boundary  at  the  edges  will  not  cause  significant  deterioration  in 
#  the  accuracy  of  the  algorithm. 


The  absorbing  boundary  conditions  are  applied  to  the  electric  fields  which  are 
tangential  to  the  outer  boundary.  Applying  center  differences  to  Equation  (2.79)  for 
the  -\-x  boundary  (i  =  /+),  and  then  using  temporal  and  spatial  averages  to  obtain 
Ey,  and  at  positions  and  times  where  they  are  calculated,  yields 


+  5. ^)  =  -  hi  +  Ik) 

J  +  5.  k)  +  +  1,  *)] 
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-2  [bJ(/+,  i  +  1.  i)  +  £,"(/+  -  1,  j  +  1,  A)] 
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2Az‘^{At  +  Ax) 
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+  £”(/+,  j  +  l.fc  -  1)  +  f^(/+  -  l,j  +'^,k-  1)} 


(2.81) 


E;^Hl*,hk+  =  -Er'(I*  -  IJ,  k  +  i) 
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(2.82) 
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When  Ey  lies  on  both  a  perfectly  conducting  ground  plane  and  the  absorbing  bound¬ 
ary,  it  must  simply  satisfy  the  boundary  condition  for  the  perfect  conductor,  and  is 
set  to  zero.  However,  when  lies  just  above  a  perfectly  conducting  ground  plane  and 
on  the  absorbing  boundary,  image  points  are  used  in  conjunction  with  the  second- 
order  absorbing  boundary  condition.  The  difference  equation  for  Ez  at  ^Az  above  a 
perfectly  conducting  ground  plane  lying  in  the  x-y  plane  at  z  —  kAz,  and  on  the  -ha: 
face  absorbing  boundary  {i  =  I'^),  is 


E;*Hl+J,k  +  i)  =  -  l,j,k+k) 

At  —  Ax 
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At  -h  Ax 
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2A!,2(A7-  +  A^)  +  5)  +  -  1.^'  +  1.  <=  +  5) 


At  +  Ax 
Ar^Ax 


-2 
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(2.83) 


(2.84) 


Applying  center  differences  to  Equation  (2.80)  at  the  +y  edge  of  the  -hx  face,  for 
Ey  ((i  =  I'^,j  =  J+  —  |)  and  Ez  {i  =  =  J'*'  —  1),  and  then  using  temporal  and 
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spatial  averages  to  obtain  Ey  and  E^  at  positions  and  times  where  they  are  calculated, 
yields 


E^+\I+,  J+  -l,k)=  E^{I+  -  1,  -  I  k)  (2.85) 

k)  +  e;{i\  -  i, k)] 

E^+\I+,  J+,  A:  +  i)  =  E^(I+  -  1,  J+,  A;  +  i)  (2.86) 

'= + ')] 

The  difference  equations  satisfying  the  ABC  for  the  tangential  electric  fields  at 
the  other  edges  and  faces  of  the  computational  domain  have  the  same  form,  and  hence 
are  not  shown  here. 

The  absorbing  boundary  is  not  required  to  absorb  the  incident  plane  wave, 
because  the  scattered  fields,  from  which  the  incident  plane  wave  has  been  analytically 
removed,  are  calculated  outside  the  inner  region.  With  infinite  ground  planes,  the 
reflected  plane  wave  is  also  analytically  removed  from  the  scattered  field  region  and 
need  not  be  absorbed  by  the  absorbing  boundary.  Since  the  scattered  field  radiates 
from  the  scatterer,  the  outgoing  scattered  fields  will  not  be  incident  on  the  absorbing 
boundary  at  grazing  angles,  which  allows  the  ABC  to  absorb  a  large  percentage  of 
the  outgoing  scattered  fields. 
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2.12  Calculation  of  Radiation  or  Scattering  Pat¬ 
terns 


The  calculation  of  the  radiation  or  scattering  patterns  requires  the  use  of  Huygens’ 
principle.  Huygens’  principle  states  that  the  field  solution  in  a  region  is  completely 
•  determined  by  the  tangential  fields  on  a  surface  which  completely  surrounds  the 

region  [106].  The  mathematical  formulation  of  Huygens’  principle  for  free  space  in 
three-dimensions,  assuming  an  time  dependence,  has  the  following  forms: 


E{f)  =  f^^dS'{iunoG{f,f')-hxH{f')  +  VxG{f,f')-nxE{f')}  (2.87) 
H{f)  =  dS'  {-iueoW,  f')  •  n  X  E{f')  +  V  x  ^(f,  f')  •  n  x  H{f')}  (2.88) 


r  =  XX  +  yy  +  zz  =  rr  (2.89) 

r'  =  xx'  +  yy'  +  iz'  =  fr'  (2.90) 


r  =  yjx^  -t-  (2-91) 

e  =  tan-i-^^=  (2.92) 

-f  y^ 

(f)  =  tan”^  —  (2.93) 

y 

f  =  £  sin  0  cos  (f)  +  y  sin9  sincf)  +  z  cos  9  (2.94) 

9  =  X  cos  9  cos  (l>  +  y  cos  9  sin  (j)  —  z  sin  9  (2.95) 
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(f)  —  —xsm(f)  +  y  cos  0 


(2.96) 
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/  +  ivv 
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giA:o|r-r'| 

J  47r  |f  —  f' 


(2.97) 


In  the  far-field,  V  can  be  approximated  as  ikof. 


_  £/■  f  r=  1  e**^o|r-r'  _ 

E{f)  =  £dS'|ic.;,.[7-rr]— 


gtfco|r-r'|  _ 
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r,,  Jko\r-f'\ 
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(2.98) 


(2.99) 


In  the  finite  difference-time  domain  scheme,  it  is  relatively  simple  to  obtain 
the  fields  over  a  closed  surface.  For  a  finite  geometry,  the  scattered  fields  can  be 
sampled  on  a  rectangular  box  which  encloses  the  entire  geometry  [Figure  2.6].  For 
an  infinite  ground  plane  geometry,  the  scattered  fields  are  sampled  on  the  sides  of  a 


2.12.  CALCULATION  OF  RADIATION  OR  SCATTERING  PATTERNS 


81 


box  which  encloses  any  discontinuities  in  or  above  the  infinite  ground  plane  [Figure 
2.7].  Image  theory  allows  the  replacement  of  the  infinite  ground  plane  with  image 
sources.  Using  the  appropriate  image  fields  for  infinite  perfect  electric  and  magnetic 
conducting  planes,  it  is  possible  to  obtain  the  tangential  fields  over  a  closed  surface. 
Note  that  the  field  solution  obtained  using  the  image  sources  is  not  valid  in  the  image 
half  space. 

Sinusoidal  and  Gaussian  pulse  time  dependent  excitations  are  treated  differently 
in  order  to  obtain  the  necessary  time  harmonic  complex  amplitudes  of  the  fields  on 
a  closed  surface.  For  sinusoidal  time  dependent  excitations,  it  is  relatively  simple 
to  obtain  the  complex  amplitudes  of  the  fields  by  sampling  the  fields  after  steady 
state  has  been  reached.  The  amplitudes  can  be  obtained  by  recording  the  maximum 
values  of  the  fields.  The  relative  phases  can  be  obtained  by  recording  and  comparing 
the  time  of  the  maximum  values  of  the  fields.  For  Gaussian  pulse  time  dependent 
excitations,  Fourier  transformation  of  the  excitation  source  and  the  fields  is  performed. 
At  a  particular  frequency,  by  dividing  the  complex  Fourier  amplitudes  of  the  fields 
by  the  Fourier  amplitude  of  the  excitation,  the  complex  amplitudes  of  the  fields 
can  be  obtained  for  an  excitation  source  with  unity  amplitude.  For  sinusoidal  time 
dependent  excitations,  fields  and  hence  radiation  patterns  can  only  be  obtained  for  a 
single  frequency.  However,  for  Gaussian  pulse  time  dependent  excitations,  fields  and 
radiation  patterns  can  be  obtained  for  multiple  frequencies.  Some  limitations  to  the 
range  of  frequencies  that  can  be  analyzed,  are  that  the  grid  must  be  fine  enough  to 
adequately  model  the  frequency  of  interest  and  that  the  Gaussian  pulse  must  contain 
a  significant  amount  of  energy  at  that  frequency. 
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2.13  Limitations 

While  the  FD-TD  technique  can  provide  very  accurate  predictions  and  solutions  to 
various  electromagnetic  phenomena,  there  are  a  number  of  approximations  inherent  to 
the  technique.  The  center  difference  approximations  applied  to  Maxwell’s  equations 
are  accurate  to  the  second  order  (i.e.,  the  error  term  is  proportional  to  the  square  of  the 
grid  size  (A^)).  This  error  can  be  kept  to  a  minimum  by  choosing  the  grid  size  to  be 
less  than  a  tenth  of  the  shortest  wavelength  of  interest.  Some  dispersion  is  introduced 
by  the  discretization  of  Maxwell’s  equations  (i.e.,  different  frequency  components 
travel  at  slightly  different  velocities).  Geometries  are  discretized  on  a  rectangular  grid. 
Hence,  the  dimensions  of  any  scatterer  are  restricted  to  integer  multiples  of  the  grid 
size,  and  curved  surfaces  must  be  approximated  with  staircases.  In  general,  as  long 
as  the  grid  size  is  small  compared  to  the  dimensions  of  the  scatterer,  the  scatterer  can 
be  adequately  modeled  on  rectangular  grids.  At  the  outer  boundary  some  reflections 
will  occur  because  the  absorbing  boundary  condition  is  an  approximate  condition. 
These  reflections  can  be  kept  tolerable  by  locating  the  outer  boundary  at  least  half  a 
wavelength  away  from  the  scatterer  [10]. 

2.14  Summary 

The  main  features  of  the  FD-TD  technique  are  the  following.  It  is  essentially  an 
“exact”  numerical  method,  with  most  of  the  errors  being  introduced  to  the  finite 
discretization  of  space  and  time.  The  method  is  well  suited  for  complex  geometries. 
Since  it  is  a  time  domain  method,  the  fields  can  be  viewed  as  a  function  of  time,  and 
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through  the  use  of  Fourier  transforms,  multiple  frequencies  can  be  analyzed  with  a 
single  FD-TD  simulation. 


CHAPTER  2. 
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Chapter  3 


EM  Radiation  from  a  VLSI  Chip 
Package  and  Heatsink  Structure 

3.1  Introduction 


The  extension  of  regulatory  electromagnetic  emissions  limits  to  the  gigahertz  range 
of  frequencies  for  high  performance  computing  equipment  has  necessitated  caution 
in  the  design  and  implementation  of  components,  down  to  the  VLSI  chip  package 
level,  for  computing  systems.  The  avoidance  of  potentially  efficient  radiators  at  the 
design  and  development  phases  is  of  interest.  The  objective  is  to  minimize  the  need 
for  cost-prohibitive  corrective  measures  commonly  invoked  in  reaction  to  problems 
encountered  in  tests  on  assembled  systems. 

In  this  chapter  the  effects  of  a  heatsink,  placed  over  a  chip  package,  on  the 
electromagnetic  emissions  characteristics  of  the  chip  and  package  configuration  are 
addressed  [31]-[33].  Such  an  interest  is  justified  by  the  following  reasons.  The  heatsink 
is  typically  metallic  and  is  in  close  proximity  to  a  major  energy  source  —  the  chip. 
The  radiation  characteristics  of  the  heatsink  will  be  similar  to  those  of  a  microstrip 
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patch  antenna  [116]-[118],[111].  The  electrical  dimensions  of  the  heatsink  structure  are 
comparable  to  the  wavelengths  of  concern.  In  light  of  the  ever-increasing  power  levels 
of  high  performance  integrated  circuitry,  the  heatsink  is  an  indispensable  component 
of  the  VLSI  chip  packaging  configuration.  The  objectives  are  to  identify  and  quantify 
the  effects  of  the  heatsink  on  the  radiating  properties  of  the  package  structure;  to 
understand  and  exploit  the  radiation  mechanisms;  and  to  evaluate  viable  heatsink 
implementation  schemes  for  minimizing  the  overall  radiating  capacity. 

In  Section  3.2,  the  model  used  to  analyze  the  problem  is  presented.  The  method 
of  solution,  namely  the  finite-difference  time-domain  technique,  is  discussed  in  Section 
3.3.  In  Section  3.4,  results  are  shown  for  two  example  heatsinks.  Observations  are 
made  regarding  the  general  radiation  characteristics  of  the  VLSI  chip  package  and 
heatsink  structure.  Two  implementation  schemes  to  reduce  the  radiation  are  also 
explored  and  evaluated.  Overall  observations  are  made  in  Section  3.5  regarding  the 
effectiveness  of  the  analysis  technique  and  the  radiation  characteristics  of  the  heatsink 
and  chip  package. 
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3.2  Problem  Configuration  and  Model 

The  electromagnetic  radiation  properties  of  a  heatsink  and  integrated  circuit  (IC) 
package  configuration  are  analyzed  by  means  of  suitable  models.  A  typical  heatsink/ 
package  configuration  of  interest  is  shown  in  Figure  3.1.  The  heatsink  may  be  mod¬ 
eled  as  a  perfectly  conducting  rectangular  slab  positioned  over  a  finite-size  dielectric 
medium  representing  the  chip  package.  The  chip  package  is  supported  in  turn  by 
a  dielectric  layer  of  infinite  extent  over  an  infinite  ground  plane  which  models  the 
substrate  or  printed  wiring  board  (PWB)  with  (at  least)  one  reference  layer.  The 
electromagnetic  sources  used  to  model  the  active  chip  include  electric  and  magnetic 
dipoles  of  vertical  and  horizontal  orientations  positioned  on  the  package/substrate 
dielectric  interface.  The  choice  of  such  elements  circumvents  the  problem  of  deriving 
a  rigorous  model  for  the  VLSI  circuitry  as  an  electromagnetic  source  and  eliminates 
factors  which  may  otherwise  clutter  the  understanding  and  interpretation  of  the  phys¬ 
ical  processes  associated  with  the  spurious  radiation  study.  The  resulting  model  of  an 
electromagnetically-coupled  slab  structure  is  shown  in  Figure  3.2.  The  use  of  a  model 
with  reduced  complexity  for  both  source  and  heatsink  is  to  facilitate  the  interpretation 
of  results  based  on  underlying  physics.  The  extension  to  other  heatsink  configurations 
(e.g.,  circular  or  finned)  and  source  models  may  be  accomplished  with  the  numerical 
solution  employed.  With  sufficiently  fine  grids,  circular  heatsinks  may  be  adequately 
modeled.  As  illustrated  in  Section  3.4,  the  model  in  Figure  3.2  is  readily  modified 
to  model  various  heatsink  implementation  schemes  proposed  for  emissions  reduction 
consideration.  Field  strengths,  and  radiated  powers  are  computed  for  quantitative 
analysis  and  evaluation. 
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Figure  3.1:  Cut  away  view  of  a  VLSI  package/heatsink  configuration. 


SIDE  VIEW 


HEATSINK  (CONDUCTING  SLAB) 


REFERENCE  PLANE 


Figure  3.2:  Cross  sectional  view  of  simplified  package/heatsink  model 
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3.3  Method  of  Solution 

The  method  employed  in  the  analysis  is  the  finite-difference  time-domain  (FD-TD) 
technique  [1]-[12].  A  more  in  depth  discussion  of  the  FD-TD  technique  is  provided  in 
Chapter  2.  It  is  based  on  the  discretization  of  the  electric  and  magnetic  fields  over 
rectangular  grids  together  with  the  finite  difference  approximation  of  the  spatial  and 
temporal  derivatives  appearing  in  the  differential  form  of  Maxwell’s  equations.  The 
reasons  for  which  the  FD-TD  methodology  was  selected  include  the  relative  ease  of 
implementation  for  complicated  geometries  with  dielectrics,  the  requirement  of  only 
simple  arithmetic  operations  in  the  solution  process,  and  the  flexibility  for  time-  and 
frequency-domain  analyses. 

In  the  FD-TD  technique,  a  computational  domain  is  first  defined  and  divided 
into  rectangular  cells.  Electric  and  magnetic  fields  are  spatially  discretized  in  a  stag¬ 
gered  manner  [1].  Electric  fields  are  assigned  to  half-integer  (n-|- 1/2)  time  steps  and 
magnetic  fields  are  assigned  to  integer  (n)  time  steps  for  the  temporal  discretization 
of  fields.  Next,  the  spatial  and  temporal  derivatives  of  the  two  Maxwell’s  curl  equa¬ 
tions  are  approximated  using  center  differences.  Maxwell’s  divergence  equations  are 
ignored  since  the  curl  equations  with  appropriate  boundary  conditions  uniquely  de¬ 
termine  the  solution.  In  rectilinear  coordinates,  the  curl  equations  can  be  rewritten  as 
a  set  of  six  scalar  equations.  Difference  equations  are  derived  from  these  six  equations 
by  applying  center  differencing.  The  center  difference  ensures  that  the  spatial  and 
temporal  discretizations  are  of  second-order  accuracy,  where  errors  are  proportional 
to  the  square  of  the  cell  size  and  time  increment  [1].  Finally,  with  appropriate  ini¬ 
tial  and  boundary  conditions,  the  solutions  to  the  difference  equations  are  obtained 
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through  explicit  leapfrog  time  marching.  This  corresponds  to  alternating  the  advance 
of  electric  and  magnetic  fields  [6].  To  achieve  accurate  results,  the  cell  sizes  are  taken 
to  be  a  fraction  (~  1/20)  of  the  smallest  wavelength.  The  time  increment  and  the 
cell  size  are  related  by  the  stability  criterion  [6], 

At  <  — ,  ^  (3.1) 

Cov^(l/Aa;)2  +  (1/Ay)2  +  {1/Az)2 

where  cq  is  the  speed  of  light  in  free  space.  Fields  are  set  to  be  zero  initially  everywhere 
to  satisfy  the  causality  condition  consistent  with  zero  excitation  for  time  less  than 
zero.  The  boundary  conditions  are  continuity  of  tangential  electric  and  magnetic  fields 
on  material  interfaces,  vanishing  tangential  electric  fields  on  perfect  conductors,  and 
the  absorbing  boundary  conditions  (ABC)  on  the  boundary  of  the  computational 
domain  [67]- [75].  First  and  second-order  absorbing  boundary  conditions  [67]  are  used 
to  limit  the  computational  domain  by  simulating  unbounded  space. 

The  size  of  the  discretization  cells  should  be  reasonably  small  in  order  to  model 
features  of  the  heatsink  and  capture  correct  field  variations.  On  the  other  hand, 
the  heatsink  should  not  be  too  close  to  the  computational  boundary  in  order  to  avoid 
significant  reflections  from  the  ABC.  The  cell  count  in  the  region  outside  the  heatsink 
will  be  very  large  when  a  small  cell  size  is  used  with  uniform  gridding.  To  avoid  the 
excessive  number  of  cells  which  would  otherwise  be  generated,  a  three-zone  gridding 
scheme  [10],[31]-[33]  is  used  (Figure  3.3).  The  first  zone  contains  the  heatsink  and 
dipole  source;  it  has  the  finest  grid  among  the  three  zones.  The  second  zone  has 
the  same  height  as  the  first  and  extends  horizontally  from  the  heatsink  to  the  outer 
boundary  of  the  computational  domain.  The  cells  in  this  zone  have  the  same  vertical 
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dimension  as  those  in  the  first  zone  in  order  to  properly  resolve  the  PWB  structure. 
However,  these  cells  have  horizontal  dimensions  which  may  be  a  few  times  larger 
than  those  of  the  cells  in  the  first  zone.  The  remainder  of  the  computational  domain 
belongs  to  the  third  zone,  where  the  cells  have  the  same  horizontal  dimensions  as  the 
cells  in  the  second  zone,  and  have  a  vertical  dimension  which  may  be  a  few  times  larger 
than  the  vertical  dimension  of  the  cells  in  the  first  or  second  zones.  Field  quantities  at 
nodes  on  an  interface  between  two  zones  are  calculated  via  a  combination  of  parabolic 
curve  fitting  and  linear  interpolation,  which  maintains  the  second-order  accuracy  of 
the  algorithm  [10], [31], [32].  For  example,  in  Figure  3.4,  the  partial  derivative  of  Hy 
with  respect  to  the  vertical  direction  z  at  position  of  Ex{i,jyk)  at  the  nth  time  step 
is  given  by 


OH, 


dz 


a  —  3 


A{l  +  a)  y 


Hy{iJ,k-  1) 


1  “  a 

^A(3  +  a) 


+  13(2^  —  9  .  2 
^24(a+l)(a  +  3)^ 


IflH, 


dz^ 


(3.2) 


where  a  is  the  ratio  between  the  cell  sizes  in  the  two  regions.  As  seen  from  the 
last  term  in  the  above  equation,  second  order  accuracy  is  maintained  everywhere. 
Although  the  discretization  scheme  is  second-order  in  all  three  zones,  the  proportion¬ 
ality  constants  for  the  discretization  errors  are  different  due  to  different  cell  sizes. 
Reasonable  accuracy  is  maintained  due  to  the  fact  that  the  spatial  variation  of  the 
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Figure  3.3:  A  three-zone  gridding  scheme  for  multi-zone  gridding.  The  heatsink  is  in 
zone  1. 

fields  outside  the  first  zone  is  gradual  and  that  all  cell  sizes  remain  a  fraction  of  the 
smallest  wavelength.  The  size  of  the  time  step  is  primarily  determined  by  the  smallest 
grid  dimension.  For  simplicity,  the  variable  time  step  implementation  [51]  is  not  used. 

The  implementation  of  the  electric  and  magnetic  dipoles  of  vertical  and  hori¬ 
zontal  orientation  in  the  FD-TD  scheme  requires  normalization.  Since  the  dipoles  are 
implemented  by  exciting  the  fields  at  a  single  node  in  the  rectangular  grid,  the  di¬ 
mensions  of  the  dipole,  and  hence  their  dipole  moments,  are  not  precisely  defined.  In 
addition  to  performing  a  FD-TD  simulation  of  the  heatsink  configuration,  a  simula¬ 
tion  is  done  using  the  same  grid  dimensions  for  the  same  dipole  in  an  infinite  parallel 
plate  structure.  The  radiated  power  for  a  dipole  in  a  parallel  plate  waveguide  can 
be  analytically  related  to  its  dipole  moment  [106].  So,  the  radiated  power  from  the 
second  FD-TD  simulation  is  used  to  obtain  the  dipole  moment,  which  in  turn  is  used 
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I  (S)  I  0  I  <S) 


Hy(iJ,k-2)  Hy(iJ,k-l)  Hy{iJ,k) 


Figure  3.4:  Node  and  field  designations  for  multi-zone  gridding  in  vertical  direction. 

to  normalize  the  dipole  in  the  first  simulation. 

In  order  to  obtain  certain  frequency  domain  parameters,  such  as  the  radiated 
power  or  the  radiation  pattern,  the  complex  field  amplitudes  on  a  closed  surface  must 
be  obtained.  The  post-processing  of  full  time  domain  information  for  the  purpose 
of  generating  frequency  domain  data  is  expensive  in  storage  and  is  inefficient.  The 
complex  field  amplitudes  are  instead  calculated  at  selected  multiple  frequencies.  To 
obtain  multi-frequency  data,  the  time  waveforms  of  the  dipole  excitations  are  chosen 
to  be  modulated  Gaussian  pulses,  i.e.,  ~  cos(u}o{t  —  to))6xp[— ((t  —  to)/T)^],  where 
uo  is  the  center  frequency  of  the  excitation,  to  is  the  delay,  and  T  determines  the 
pulse-width.  The  Gaussian  pulse  is  modulated  in  order  to  concentrate  the  energy  in 
the  frequency  range  of  interest.  Complex  amplitudes  at  each  frequency  are  calculated 
simultaneously  on  a  selected  surface  outside  the  first  zone  using  the  discrete  Fourier 
transform.  The  total  radiated  power  is  found  by  integrating  the  Poynting  vector  over 
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a  closed  surface.  Once  the  complex  field  amplitudes  on  a  closed  surface  have  been 
stored,  the  far-field  radiation  pattern  can  be  calculated  using  Huygens’  principle  [106]. 
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3.4  Results  and  Discussion 

The  basic  model  described  in  Figure  3.2  is  analyzed  with  the  use  of  horizontal  electric 
and  magnetic  dipoles  and  vertical  electric  and  magnetic  dipoles  as  sources.  Of  these 
sources  only  the  horizontal  magnetic  dipole  (HMD)  and  vertical  electric  dipole  (VED) 
show  radiation  enhancement  in  the  presence  of  the  heatsink.  The  frequency  range 
considered  does  not  exceed  10  GHz.  The  dominant  horizontal  (electric  field)  polar¬ 
ization  generated  by  the  other  2  sources  remains  below  cutoff  within  this  frequency 
range  and  for  the  typical  dimensions  considered.  Both  classes  of  sources  are  found 
in  VLSI  chip  package  configurations  —  the  VED  in  the  form  of  current  carrying  vias 
and  pins  and  the  HMD  as  vertically  configured  signal  and  return  paths.  Results  and 
discussions  will  focus  on  the  VED  and  HMD  sources. 

As  stated  previously,  the  choice  of  the  dipoles  for  sources  facilitates  the  under¬ 
standing  of  the  overall  problem  and  the  interpretation  of  the  simulation  results  while 
deferring  the  electromagnetic  modeling  of  a  VLSI  chip.  From  a  practical  viewpoint, 
the  drawback  of  such  models  lies  in  the  resulting  difficulty  for  designers  to  determine 
the  radiation  levels  of  the  actual  heatsink  and  package  configurations.  A  knowledge  of 
the  electromagnetic  source  characteristics  has  to  be  combined  with  a  physical  model 
(e.g.,  a  conducting  patch)  of  the  VLSI  chip  to  meet  this  practical  requirement. 

3.4.1  Radiation  Enhancement  by  Heatsinks 

The  configuration  of  heatsink  A  is  the  same  as  the  heatsink  illustrated  in  Figure  3.2, 
except  that  no  dielectrics  are  present.  The  dimensions  of  the  square  conducting  slab 
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are  Whs  =  4.8  cm,  4*  =  2.5  mm,  and  h  =  7.5  mm,  where  Whs  is  the  width,  4^  is  the 
thickness,  and  h  is  the  height  of  the  heatsink  above  the  ground  plane.  The  dipole 
source  is  located  within  the  three-dimensional  rectangular  grid  at  the  node  closest  to 
the  center  of  the  heatsink/ground  plane  cavity. 

Figure  3.5  compares  the  normalized  radiated  power  of  an  HMD  over  a  ground 
plane  in  the  presence  and  absence  of  heatsink  A  for  a  frequency  range  of  1  to  6.5 
GHz.  Normalization  is  with  respect  to  the  total  radiated  power  of  an  identical  dipole 
in  unbounded  space  (no  heatsink  or  ground  plane).  The  resonant  behavior  around 
2.3  and  5.5  GHz  in  the  presence  of  the  heatsink  is  evident,  with  the  first  resonance 
showing  significant  enhancement  of  radiation  over  the  case  without  the  heatsink.  The 
3  dB  level  plotted  corresponds  approximately  to  the  normalized  dipole  radiation  with 
the  ground  plane  but  in  the  absence  of  the  heatsink. 

It  may  be  inferred  from  the  heatsink  dimensions  and  the  field  patterns  gen¬ 
erated  by  a  centered  HMD,  that  the  resonant  frequency  at  approximately  2.3  GHz 
corresponds  to  the  lowest  resonance.  The  electric  fields  generated  by  a  HMD  are 
anti-symmetric  in  directions  orthogonal  to  the  dipole  axis.  Hence,  a  centered  HMD 
will  only  excite  modes  with  at  least  one  degree  of  anti-symmetry  in  the  horizontal 
direction  which  is  orthogonal  to  the  dipole  axis. 

Figure  3.6  shows  the  corresponding  results  over  the  frequency  range  of  1  to  8 
GHz  with  a  centered  VED  as  the  source.  Again,  substantial  enhancement  in  radiated 
power  in  the  presence  of  a  heatsink  occurs  at  the  primary  resonance  of  4.9  GHz,  and 
at  the  second  resonance  near  7  GHz.  The  electric  fields  generated  by  the  VED  are 
symmetric  about  the  dipole  axis,  hence  only  symmetric  cavity  modes  are  expected  to 
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NORMALIZED  RADIATED  POWER  WITH/WITHOUT  HEATSINK 

NORMALIZED  TO  DIPOLE  RADIATION 

Figure  3.5:  Comparison  of  normalized  total  radiated  power  in  presence  and  absence 
of  heatsink  A  with  a  HMD  source. 
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NORMALIZED  RADIATED  POWER  WITH/WITHOUT  HEATSINK 

NORMALIZED  TO  DIPOLE  RADIATION 


Figure  3.6:  Comparison  of  normalized  total  radiated  power  in  presence  and  absence 
of  heatsink  A  with  an  VED  source. 


be  excited. 

The  enhancement  of  the  total  radiated  power  over  when  the  heatsink  is  absent 
in  the  vicinity  of  the  resonance  is  evident.  In  particular,  for  the  lowest  resonant 
frequency  near  2.3  GHz,  16  dB  enhancement  is  observed  with  the  addition  of  the 
heatsink  to  the  ground  plane.  The  actual  degree  of  enhancement  will  be  influenced 
by  the  physical  chip  model,  the  module  and  system  environment.  It  is  observed,  in 
the  cases  shown,  that  the  resonances  are  not  sharp,  i.e.,  they  are  of  relatively  low 
quality  factor,  Q.  For  instance,  in  Figure  3.5,  the  radiation  enhancement  is  higher 
than  10  dB  over  a  frequency  range  of  2.05  GHz  to  2.55  GHz.  The  500  MHz-wide  band 
could  fit  4-5  harmonics  of  a  100  MHz  clock  signal,  thereby  enhancing  radiation  at  a 


3.4.  RESULTS  AND  DISCUSSION 


99 


multitude  of  frequencies  in  general.  An  approximate  determination  of  the  resonant 
frequencies  is  of  interest  and  warrants  further  discussion. 

A  simple  but  less  accurate  method  of  determining  the  resonant  frequencies  for 
such  a  configuration  is  the  magnetic-wall  approximation  [116], [111].  The  associated 
formula  is  determined  from  the  assumption  of  perfect  open-circuited  boundaries  and 
also  assumes  very  small  cavity  height.  The  formula  for  the  frequency  associated  with 
the  (m,  n)  mode  takes  the  form 


where  c  is  the  speed  of  light  within  the  cavity  medium.  In  the  above  equation,  m  and 
n  are  used  to  denote  the  mode  index  corresponding  to  the  m  (n)  half-wavelength  field 
variation  in  the  x  (y)  direction.  Because  the  cavity  height  is  small  at  all  frequencies 
of  interest,  the  third  index  for  field  variation  in  the  vertical  direction  is  not  used, 
and  is  implicitly  zero.  and  Wy  are  the  effective  horizontal  dimensions  derived 
from  enlarging  the  width  and  depth  to  reflect  the  effect  of  fringing  fields.  Such 
enlargement  formulas  abound  in  microstrip  literature  [111],  and  are  dependent  on 
conductor  thickness,  cavity  height  and  the  dielectric  environment.  No  simple  formulas 
are  available  which  include  the  effects  of  the  corners. 

Table  3.1  lists  the  first  few  resonances  for  the  heatsink  configuration  examined 
in  Figures  3.5  and  3.6  as  predicted  by  the  magnetic-wall  formula  using  effective  di¬ 
mensions.  Using  the  dimensions  for  heatsink  A,  the  effective  width  {Wx)  and  depth 
(kFj,)  are  approximately  5.0  cm.  Also  shown  are  the  corresponding  resonant  frequen- 


100 


CHAPTER  3.  RADIATION  FROM  HEATSINK  STRUCTURE 


Mode 

m  n 

Resonant  Frequency  fm,n 
Magnetic- Wall  (GHz)  FD-TD  (GHz) 

FD-TD  Source  Type 

0 

1 

3.0 

2.3 

centered  HMD 

0 

2 

6.0 

4.9 

centered  VED 

1 

1 

4.3 

3.4 

off-center  HMD 

1 

2 

6.7 

5.5 

centered  HMD 

2 

2 

8.5 

7.0 

centered  VED 

Table  3.1;  Magnetic-wall  resonances  and  corresponding  FD-TD  predicted  resonant 
frequencies  for  heatsink  A. 

cies  derived  from  the  source-dependent  FD-TD  simulations.  In  particular,  the  first 
resonant  mode  corresponds  to  the  longer  effective  dimension  being  equal  to  one-half 
of  the  wavelength  of  the  associated  resonant  frequency  with  no  field  variation  in  the 
other  dimension.  This  may  be  excited  by  a  HMD  located  at  the  center  of  the  cavity. 
The  lowest  resonance  excited  by  a  centered  VED  is  the  (2,0)  mode  which  corresponds 
to  the  longer  dimension  having  1-wavelength  variation  and  no  variation  in  the  other 
dimension. 

Noted  are  the  discrepancies  in  values  attributed  to  the  low  height  and  limited 
fringing  field  constraints  on  the  magnetic- wall  solution.  Also  notable  is  the  absence 
of  correspondence  in  the  FD-TD  simulations  illustrated  in  Figures  3.5  and  3.6  for  the 
(1,1)  magnetic-wall  resonance.  This  resonance  is  excited  by  off-center  dipoles  and 
will  be  illustrated  in  Figure  3.8. 

For  the  heatsink  configuration  examined  in  Figure  3.5,  Figure  3.7  shows  the 
surface  plots  of  the  vertical  electric  field  amplitudes  at  4  different  time  steps  over 
a  horizontal  plane  that  traverses  the  heatsink  and  ground  plane  cavity  and  spans 
the  computational  domain.  Each  plot  is  normalized  to  the  maximum  value  for  the 
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respective  time  step.  The  absence  of  any  obvious  reflected  disturbances  from  the  edges 
of  the  computational  domain  illustrates  the  efficiency  of  the  absorbing  boundary 
conditions  implemented.  The  sustained  resonant  mode  pattern  corresponds  to  the 
(1,0)  mode  and  is  evident  after  1000  time  steps. 

Heatsink  B,  which  is  a  little  larger  than  heatsink  A  is  also  analyzed.  Heatsink 
B  has  dimensions  of  Wh$  =  5.8  cm,  4s  =  1-0  cm,  and  h  =  1.3  cm  (Figure  3.2)  and  this 
model  contains  no  dielectrics.  The  dipole  source  is  again  located  at  the  node  closest 
to  the  center  of  the  heatsink/ground  plane  cavity  unless  stated  otherwise. 

Figure  3.8  compares  the  normalized  total  radiated  power  of  an  HMD  over  a 
ground  plane  in  the  presence  and  absence  of  heatsink  B  for  a  frequency  range  of  1  to 
3  GHz.  The  first  resonance  occurs  at  1.75  GHz,  and  shows  radiation  enhancement  of 
approximately  13  dB  over  when  the  heatsink  is  absent.  In  comparison  to  the  normal¬ 
ized  radiated  power  with  heatsink  A,  the  normalized  radiated  power  with  heatsink 
B  has  a  lower  resonant  frequency  and  has  a  broader  response  in  frequency.  This  is 
consistent  with  the  larger  dimensions  of  heatsink  B  and  its  larger  displacement  from 
the  ground  plane  compared  to  heatsink  A  [116], [111].  Figure  3.8  also  shows  the  corre¬ 
sponding  response  for  an  off-center  HMD  source  for  the  same  physical  configuration. 
A  weaker  excitation  of  the  lowest  resonance  and  the  excitation  of  the  next  resonance 
at  2.7  GHz,  absent  for  centered  dipole  excitation,  are  manifested. 

In  considering  the  radiation  enhancement  as  a  function  of  frequency,  the  concern 
begins  with  the  lowest  resonance  which  is  excited  by  a  centered  HMD  source.  The 
remaining  illustrations  will  be  based  on  the  centered  HMD  as  the  source. 


102 


CHAPTER  3.  RADIATION  FROM  HEATSINK  STRUCTURE 


Figure  3.7:  Vertical  electric  field  amplitude  plots  over  horizontal  plane  through 
heatsink/  ground  plane  cavity  for  HMD  source.  Source  function  ~  cos(a;o(t  — 
to))exp[-((t  -  to)/T)%  where  a;o  =  Stt  x  10®  rad/s,  to  =  4T  =  0.8  ns.  N  is  the 
number  of  time  steps  elapsed,  and  each  time  step  corresponds  approximately  to  2.2 
picoseconds. 
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NORMALIZED  RADIATED  POWER  WITH/WITHOUT  HEATSINK 

CENTERED  (0.0.5)  AND  OFF-CENTER  (9.9.5)  HMD 

Figure  3.8:  Comparison  of  normalized  total  radiated  power  in  presence  and  absence 
of  heatsink  B  with  a  centered  and  an  off-center  HMD  source. 
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NORMALIZED  RADIATED  POWER  WITH  HEATSINK 

HMD  -  HEATSINK  (5.8  X  5,8  x  1.0  cm)  CAVITY  WITH  VARIABLE  HEIGHT 

Figure  3.9:  Comparison  of  normalized  total  radiated  power  with  HMD  source  for 
heatsink  5  as  a  function  of  cavity  height. 

In  Figure  3.9,  the  variation  in  the  normalized  total  radiated  power  for  heatsink 
B  with  changes  in  the  cavity  height  is  shown.  The  three  cavity  heights  considered  are 
0.6  cm,  0.9  cm  and  1.3  cm.  The  downward  shift  in  resonant  frequency  with  increasing 
cavity  height  is  consistent  with  the  increased  fringing  fields  which  results  in  larger 
effective  dimensions.  The  broadening  of  the  peaks  in  the  curves  indicates  decreasing 
Q’s  which  is  expected  for  increasing  cavity  heights. 

The  consideration  of  a  heatsink  without  dielectrics  so  far  has  allowed  the  verifi¬ 
cation  of  expected  behavior  and  trends.  It  is  appropriate  to  verify  the  performance  of 
the  algorithms  for  the  more  realistic  case  where  dielectrics  such  as  the  PWB  and  the 
chip  package  are  included.  The  dimensions  of  this  model  (Figure  3.2)  are  Whs  =  5.8 
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cm,  ths  —  1-0  cm,  h  =  1.3  cm,  tp^ub  =  2.5  mm,  and  Wic  =  5.1  cm,  Uc  =  6.4  mm. 
The  PWB  is  modeled  as  a  lossless  dielectric  with  dielectric  constant  of  4.0  and  the 
chip  package  as  a  dielectric  with  dielectric  constant  of  9.0  and  a  small  conductivity 
of  0.0001  mhos/m.  The  normalized  total  radiated  power  results  are  shown  in  Figure 
3.10.  The  expected  trend  of  lower  resonant  frequencies  and  higher  Q  in  the  presence 
of  dielectric  material  is  manifested.  The  inclusion  of  the  dielectric  material  results 
in  an  effective  dielectric  constant  larger  than  1,  leading  to  a  shorter  wavelength  than 
in  free  space.  Hence,  the  heatsink  appears  larger  electrically,  and  has  lower  reso¬ 
nant  frequencies.  It  is  noted  that  no  alternative  methods  to  the  FD-TD  method  are 
readily  available  for  the  analysis  of  the  3-dimensional  mixed-dielectric  problem  that 
results  when  the  PWB  and  the  chip  package  models  are  considered  together  with  the 
heatsink. 

A  typical  heatsink  with  fins  has  been  modeled  in  the  preceding  sections  as 
a  perfectly  conducting  rectangular  slab.  The  validity  of  such  an  approximation  is 
examined  in  the  following  discussion.  The  dimensions  of  the  finned  heatsink  are  such 
that  it  would  fit  exactly  within  heatsink  B.  The  fins  are  7.5  mm  in  height  and  are 
evenly  distributed  across  the  surface  of  the  heatsink.  The  fins  begin  at  a  distance  of 
2.5  mm  from  the  lower  surface  of  the  heatsink,  and  are  of  zero  thickness  for  simplicity 
in  modeling.  The  horizontal  dimensions  and  the  total  vertical  dimension  of  the  finned 
heatsink  are  the  same  as  the  corresponding  dimensions  of  heatsink  B. 

Figure  3.11  compares  the  normalized  total  radiated  power  for  a  centered  x- 
directed  HMD  in  the  presence  of  plain  and  finned  heatsink  models,  the  latter  of  which 
bear  fins  with  either  x  ox  y  orientations.  Emphasis  is  on  the  characteristics  near  their 
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NORMALIZED  RADIATED  POWER  WITH/WITHOUT  HEATSINK 

VARIATION  WITH  PACKAGE  AND  PWB  DIELECTRIC  MODELS 

Figure  3.10:  Comparison  of  normalized  total  radiated  power  with  HMD  source  for 
heatsink  B  only,  with  PWB  model,  and  with  PWB  and  IC  package  models. 
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RADIATED  POWER  NORMALIZED  TO  DIPOLE  RADIATION 

COMPARISON  OF  HEATSINK  MODELS  WITH  AND  WITHOUT  FINS 

Figure  3.11:  Comparison  of  normalized  total  radiated  power  with  a;-directed  HMD 
source  for  plain  and  finned  versions  of  heatsink  B. 

respective  first  resonances.  The  minor  differences  in  the  resonance  position  among 
the  3  models  are  attributable  to  the  difference  in  heatsink  features  in  the  vicinity  of 
the  dominant  fringing  fields  (at  the  ±y  faces)  corresponding  to  the  lowest  resonance. 
In  the  case  of  a  plain  heatsink  model,  the  fringing  fields  result  from  a  uniformly  thick 
conductor  (lowest  resonant  frequency  of  1.75  GHz).  For  the  model  with  fins  parallel 
to  the  y  direction,  the  dominant  fringing  fields  correspond  to  a  thinner  heatsink, 
leading  to  an  effectively  smaller  cavity  with  a  higher  resonant  frequency  (1.78  GHz). 
For  the  model  with  fins  parallel  to  the  x  direction,  the  height  of  the  fins  together  with 
the  fields  existing  in  the  finned  environment  together  results  in  an  effectively  larger 
cavity  leading  to  a  reduction  in  resonance  frequency  (1.70  GHz). 
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The  differences,  nevertheless,  appear  insignificant  and  help  support  the  appli¬ 
cation  of  the  plain  model  in  the  analyses.  As  further  demonstration  of  the  validity  of 
the  plain  model  and  to  demonstrate  the  capacity  of  the  FD-TD  method  in  facilitating 
post-processing  for  radiation  characteristics,  Figure  3.12  shows  the  gain  patterns  at 
the  respective  first  resonances  for  the  finned  and  plain  heatsink  models  excited  by  a 
centered  a;-directed  HMD.  In  addition.  Figure  3.13  displays  the  electric  field  pattern 
at  resonance  over  a  vertical  plane  for  a  dipole  moment  of  10“®  A-m,  and  at  a  distance 
of  3  meters.  The  results  for  the  HMD  radiating  in  the  absence  of  the  heatsink  are 
also  included  in  both  figures  for  reference.  The  Huygens’  surface  utilized  to  calculate 
the  gain  and  field  patterns  is  a  box  which  entirely  encloses  the  heatsink  structure. 
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Figure  3.12:  Gain  patterns  at  the  respective  resonances  (1.75,  1.70,  1.78  GHz)  of 
unfinned,  x-  and  y-directed  finned  heatsink  models  excited  by  a  centered  rr-directed 
HMD.  The  gain  pattern  for  a  HMD  (1.75  GHz)  radiation  without  a  heatsink  is  also 
included. 
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Figure  3.13:  Electric  field  patterns  at  the  respective  resonances  (1.75, 1.70,  1.78  GHz) 
of  unfinned,  x-  and  y-directed  finned  heatsink  models  excited  by  a  centered  a:-directed 
HMD.  The  distance  is  3  meters  and  the  dipole  moment  is  10“®  A-m. 
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3.4.2  Heatsink  Implementation  Options  for  Reduction  of 
Radiation 

From  the  simplified  model  for  the  heatsink  configuration  it  is  evident  that  cavity 
resonances  exist,  which  may  enhance  radiation  over  when  the  heatsink  is  absent. 
The  major  assumptions  have  been  the  neglect  of  any  adjacent  modules,  and  other 
components  on  a  module;  the  neglect  of  metalization  in  the  heatsink  cavity  attributed 
to  interconnects,  e.g.,  package  references,  leads,  and  pins;  and  the  assumption  that 
basic  dipole  sources  are  adequate  source  models. 

Under  such  circumstances,  and  given  that  the  heatsink  will  likely  continue  as 
an  indispensable  component  of  the  VLSI  packaging  environment,  it  is  of  interest  to 
explore  the  feasibility  of  some  implementation  schemes  which  exploit  the  heatsink 
presence  to  reduce  the  level  of  radiation. 

To  this  end,  two  implementation  options  are  analyzed:  the  use  of  metallic  posts 
connecting  the  edges  of  the  heatsink  to  the  reference  plane  and  the  application  of 
a  conducting  gasket  around  the  lip  of  the  heatsink/reference  plane  cavity.  Both 
configurations  are  shown  in  Figures  3.14  and  3.15. 

In  the  posting  option,  two  cases  are  considered:  the  use  of  a  single  post  at  each 
of  the  4  corners  (4-post  option)  or  at  the  corners  and  centers  of  the  edges  (8-post 
option).  Numerical,  FD-TD,  simulations  are  performed  assuming  perfect  connections 
between  the  posts,  heatsink  and  reference  plane.  In  the  gasketing  option,  the  ring 
of  gasket  material  is  excised  at  the  base  of  each  of  the  4  faces  to  model  openings 
which  may  be  required  for  interconnect  passage.  The  gasket  material  is  modeled  as 
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Figure  3.14:  Model  for  heatsink  B  with  4  support  posts  connected  to  reference  plane. 


Figure  3.15:  Model  for  heatsink  B  with  a  conducting  gasket. 
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RADIATED  POWER  NORMALIZED  TO  DIPOLE  RADIATION 

MULTIPLE  POSTS  OPTION 

Figure  3.16:  Comparison  of  normalized  total  radiated  power  with  HMD  source  for 
heatsink  B  with  4  and  8  posts  and  without  posts. 

a  conducting  dielectric. 

Detailed  dimensions  of  the  4-post  model  with  heatsink  B  are  shown  in  Figure 
3.14.  Figure  3.16  shows  the  results  from  simulations  on  the  4-  and  8-post  models. 
The  normalized  total  radiated  power  is  plotted  as  a  function  of  frequency  over  the 
range  of  1  to  3  GHz.  The  suppression  of  the  original  resonance  at  1.75  GHz  is  rather 
dramatic  with  4  posts  and  even  more  so  with  8  posts.  However,  the  suppression  is  not 
uniform  with  frequency  and  what  is  clear  from  the  results  is  that  new  resonances  have 
appeared  with  posting.  These  occur  at  higher  frequencies  —  2.6  GHz  and  greater 
than  3  GHz  respectively  for  4  and  8  posts.  This  is  of  significance  if  the  new  resonances 
may  be  shifted  beyond  the  frequency  range  where  emissions  limits  apply. 
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The  gasketed  heatsink  model  for  heatsink  B  is  shown  in  Figure  3.15.  Figure 
3.17  presents  the  results  from  simulations  on  this  configuration.  The  normalized  total 
radiated  power  is  shown  over  a  frequency  range  from  1  to  3  GHz.  The  gasketing  option 
appears  more  effective  at  higher  frequencies  and  may  have  a  slight  adverse  effect  at  the 
lower  frequency  end.  The  improvement  at  high  frequencies  is  not  unexpected  since, 
with  the  finite  material  conductivity,  the  skin  depth  decreases  with  frequency  and 
confines  as  well  as  dissipates  the  electromagnetic  energy  within  the  heatsink/reference 
plane  cavity  more  efficiently.  A  possible  explanation  for  the  low  frequency  behavior 
is  the  assumed  low  conductivity  of  the  gasket  material  leading  to  significant  current 
flow  on  the  outer  surface  of  the  gasket  material.  There  is  then  an  increased  area 
over  which  current  flows  and  contributes  to  the  radiation.  The  effect  of  varying  the 
conductivity  of  the  gasket  is  also  examined  in  Figure  3.17.  The  gradual  reduction  of 
the  lower  frequency  peak  as  gasket  conductivity  increases  is  noted.  This  is  consistent 
with  the  earlier  discussion  regarding  the  behavior  of  gasketing  at  low  frequencies. 
The  performance  at  a  lower  conductivity  of  5  mhos/m  near  1  GHz  worsens,  while 
at  a  higher  conductivity  of  60  mhos/m,  suppression  near  1  GHz  is  manifested.  For 
reference,  the  skin  depths  at  1  GHz  for  conductivities  of  5,  10,  and  60  mhos/m 
are  approximately  7,  5  and  2  mm  respectively  while  the  gasket  thickness  for  this 
configuration  is  approximately  6.5  mm.  The  skin  depth  rule  may  also  be  misleading 
since  it  is  derived  from  plane  wave  (far-field)  arguments  and  need  not  apply  in  this 
environment  where  the  gasket  material  is  so  close  to  the  energy  source. 

For  all  three  conductivities,  the  effect  of  cavity  resonance  near  1.75  GHz  is  still 
evident.  The  lossy,  dielectric,  gasket  model  perturbs  the  field  patterns  and  affects  the 
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RADIATED  POWER  NORMALIZED  TO  DIPOLE  RADIATION 

VARIATION  OF  GASKET  CONDUCTIVITY 


Figure  3.17:  Comparison  of  normalized  total  radiated  power  with  HMD  source 
heatsink  B  with  gasketing  of  varying  conductivity. 
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effective  cavity  dimensions  and  hence  the  position  of  the  resonant  frequencies. 

In  practice,  compliant  gasket  materials  suitable  for  such  applications  have  bulk 
conductivities  of  the  order  of  10^-10^  mhos/m.  For  such  conductivities,  the  skin  depth 
would  be  of  the  order  of  0.1  mm  at  5  GHz.  To  resolve  the  field  variations  across  the 
gasket  material  (1-2  grid  cells  per  skin  depth)  would  require  such  a  fine  discretization 
that  the  computation  requirements  for  the  3-D  problem  would  be  impractical  even 
with  the  non-uniform  discretization  scheme  used.  This  is  therefore  a  limitation  of  the 
solution  method. 

Since  the  posting  option  pushes  the  resonances  upward  in  frequency  where  gas¬ 
keting  functions  best,  it  is  logical  to  expect  radiation  reduction  over  a  broad  frequency 
band  when  both  posts  and  gasket  are  implemented.  Figure  3.18  illustrates  the  nor¬ 
malized  total  radiated  power  for  the  corresponding  model.  Although  the  level  of 
reduction  in  radiation  at  a  given  frequency  is  not  the  lowest  when  compared  to  the 
two  implementation  options  considered  separately,  it  is  observed  that  peak  radiation 
has  been  reduced  to  below  the  peaks  obtained  using  either  of  the  other  two  imple¬ 
mentation  options.  Total  radiation  is  also  always  lower  when  compared  to  the  case  of 
an  isolated  or  floating  heatsink  over  a  reference  plane.  The  observed  results  are  de¬ 
sirable  since  any  suppression  scheme  to  satisfy  regulatory  emissions  limits  should  be 
effective  over  a  broad  range  of  frequencies  as  excessive  radiation  at  a  single  frequency 
is  sufficient  to  result  in  test  failure. 

Even  though  the  analyzed  configurations  exhibit  resonances,  they  are  of  suffi¬ 
ciently  low  Q’s  that  fewer  than  8000  time  steps  are  necessary  for  the  analysis  of  any 
of  the  configurations  discussed.  For  configurations  having  higher  Q's,  the  required 
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NORMALIZED  TOTAL  RADIATED  POWER 

GASKET-AND-POSTS  OPTION 

Figure  3.18:  Comparison  of  normalized  total  radiated  power  with  HMD  source  for 
heatsink  B,  with  gasketing,  with  4  ground  posts,  and  with  combined  gasketing  and  4 
ground  posts. 
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time  response  can  be  obtained  by  using  an  extrapolation  scheme  based  on  Prony’s 
method  [114].  A  typical  computational  domain  may  contain  up  to  50  x  50  x  40  grid 
cells,  depending  on  the  frequency  range  of  interest  and  the  expansion  factors  selected 
for  the  non-uniform  gridding  implementation. 

The  numerical  simulations  performed  for  the  presented  data  on  radiated  powers 
take  on  the  order  of  5  to  30  CPU  hours  on  mini-computer  and  RISC-based  workstation 
systems,  including  pre-  and  post-processing.  A  significant  fraction  of  the  computation 
time  is  attributed  to  processing  the  complex  fields  at  multiple  frequencies  —  necessary 
for  frequency  domain  analyses.  The  actual  computation  times  also  depend  on  the 
available  memory  and  other  system  parameters. 
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3.5  Summary 

The  FD-TD  technique  has  been  applied  to  analyze  radiation  from  a  VLSI  heatsink 
and  package  configuration.  The  technique  was  chosen  for  its  flexibility  in  treating  this 
class  of  problems  which  include  potentially  complex  configurations  for  which  neither 
approximate  nor  analytical  methods  are  practical.  The  multi-zone  gridding  scheme 
implemented  allowed  a  high  grid  resolution  in  the  vicinity  of  the  heatsink  without 
sacrificing  overall  cell  number  thereby  allowing  increased  numerical  efficiency  for  low 
frequencies.  Additional  advantages  of  this  discretization  scheme  include  improved 
geometrical  and  material  modeling. 

The  analysis  of  the  heatsink,  a  critical  component  within  the  computer  packag¬ 
ing  environment,  illustrated  the  significance  of  resonance  due  to  appreciable  electrical 
dimensions  on  the  spurious  electromagnetic  radiation  from  the  VLSI  packaging  en¬ 
vironment.  Simulations  with  typical  heatsinks  dimensions  showed  the  occurrence  of 
resonance  in  the  low  gigahertz  range.  The  effects  of  the  presence  of  a  heatsink  on  the 
radiation  properties  of  dipole  models  have  been  explained  and  the  features  of  reso¬ 
nant  behavior  presented.  The  effectiveness  and  consequences  of  exploiting  practical 
heatsink  implementation  options  such  as  grounding  and  shielding  to  reduce  electro¬ 
magnetic  emissions  have  been  discussed  and  quantified  using  the  FD-TD  numerical 
technique. 
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Chapter  4 


EM  Radiation  from 

Modules-on-Backplane 

Configurations 


4.1  Introduction 


In  this  chapter,  the  electromagnetic  emissions  characteristics  of  a  modules-on-back- 
plane  configuration  found  in  computers  are  examined  [34].  Primary  motivation  for 
studying  the  effect  of  various  structures  present  within  a  computer  on  electromagnetic 
interference  (EMI)  is  that  computers  are  now  operating  at  faster  speeds  and  con¬ 
suming  more  power  which  results  in  significantly  increased  levels  of  power  at  higher 
frequencies.  The  modules-on-backplane  is  of  interest  because  it  is  a  common  pack¬ 
aging  hierarchy.  The  modules  are  the  popular  plug-in  printed  circuit  boards  which 
are  mounted  on  the  backplane  (motherboard).  Also  a  typical  modules-on-backplane 
configuration  has  reference  plane  dimensions  which  are  appreciable  electrically  at  fre¬ 
quencies  in  the  hundreds  of  megahertz,  and  the  structure,  as  a  whole,  may  function 
as  an  antenna.  The  objective  is  to  identify  and  quantify  the  radiation  mechanisms 
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in  a  modules-on-backplane  configuration  and  to  validate  the  numerical  model  with 
measurements. 

Section  4.2  presents  the  model  used  to  analyze  the  modules-on-backplane  con¬ 
figuration.  The  method  of  solution  is  the  finite-difference  time-domain  (FD-TD)  tech¬ 
nique,  which  is  discussed  in  Section  4.3.  In  addition,  Prony’s  method  is  outlined  in 
Section  4.3.  Section  4.4  describes  the  measurement  procedure  and  provides  some  val¬ 
idation  of  the  model  and  measurement  procedure.  Measurements  and  numerical  pre¬ 
dictions  are  presented  and  discussed  for  various  configurations  in  Section  4.5.  Overall 
observations  are  drawn  in  Section  4.6  regarding  the  effectiveness  of  the  method,  and 
the  electromagnetic  emission  characteristics  of  the  modules-on-backplane  configura¬ 
tion. 


4.2  Problem  Configuration  and  Model 

A  typical  configuration  for  the  modules-on-backplane  is  shown  in  Figure  4.1.  The 
elements  of  the  packaging  hierarchy  of  interest  in  this  study  are  the  backplane,  the 
connector  and  the  module.  The  modules  and  backplane  are  printed  wiring  boards 
with  at  least  one  reference  layer.  The  connector  contains  numerous  wires  for  carrying 
signals,  for  providing  power  to  the  various  power  planes,  and  for  grounding  the  ref¬ 
erence  planes.  Due  to  the  complexity  of  the  modules-on-backplane  configuration,  it 
is  necessary  to  simplify  the  problem  to  some  reasonable  level  of  detail.  The  modules 
and  backplanes  are  populated  with  numerous  chips  and  contain  many  signal  paths. 
However,  the  single  largest  conductor  or  signal  path  in  them  is  typically  one  of  the 
ground  or  power  planes.  These  planes,  because  of  their  larger  dimensions,  tend  to 
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enhance  radiation  at  lower  frequencies  more  than  the  other  elements  associated  with 
the  modules/backplane.  For  this  reason,  the  modules  and  backplane  will  be  mod¬ 
eled  as  perfectly  conducting  plates.  The  source  considered  is  the  transient,  inductive 
voltage  developed  across  the  connector  between  the  module  and  the  backplane.  The 
implementation  of  this  source  was  done  by  connecting  the  outer  conductor  of  a  coax¬ 
ial  line  to  the  backplane  and  the  center  conductor  to  the  module.  The  choice  of 
a  coaxial  probe  to  excite  the  module  was  made  primarily  because  of  its  simplicity. 
In  addition,  the  treatment  of  multi-pin  connectors  that  may  provide  simultaneous 
non-uniform  excitation  should  in  principle  be  analyzable  through  the  use  of  superpo¬ 
sition.  The  resulting  simplified  model  of  the  modules-on-backplane  configuration  is 
shown  m  Figure  4.2.  The  use  of  a  model  with  reduced  complexity  for  both  source  and 
modules-on-backplane  configuration  is  to  facilitate  the  interpretation  of  results  based 
on  underlying  physics.  The  extension  to  other  modules-on-backplane  configurations 
(e.g.,  multiple  modules,  enclosed)  may  be  accomplished  with  the  numerical  method 
employed.  In  addition,  as  other  significant  excitation  sources  are  identified  and  better 
understood,  they  can  be  incorporated  into  the  simulation  model. 
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4.3  Method  of  Solution 

4.3.1  FD-TD  Technique 

The  finite  difference-time  domain  (FD-TD)  technique  is  used  to  analyze  this  problem 
[1]-[12].  This  method  was  chosen  because  it  has  the  capability  of  analyzing  multiple 
frequencies  with  a  single  simulation,  is  well  suited  for  low  loss  models  and  rectilinear 
geometries,  and  is  a  proven  method,  as  it  has  been  previously  used  to  successfully 
analyze  radiation  from  a  VLSI  chip  package  and  heatsink  structure  [31]-[33]. 

A  rectangular  three-dimensional  grid  is  used  as  the  computational  domain.  Spa¬ 
tially,  the  computational  domain  is  discretized  such  that  the  largest  grid  dimension 
is  less  than  about  a  twentieth  of  the  shortest  wavelength  of  interest.  A  higher  res¬ 
olution  was  chosen  over  the  typical  tenth  of  a  wavelength  criteria  in  order  to  more 
accurately  model  the  resonance  phenomena  associated  with  this  configuration.  The 
time  step  approximately  must  be  less  than  the  smallest  grid  dimension  divided  by 
the  speed  of  light  in  order  to  guarantee  the  stability  of  the  algorithm  [6].  Some 
cell  dimensions  which  were  typical  of  the  simulations  performed  for  this  paper  are 
Ax  =  5.00  mm,  Ay,=  3.65  mm,  and  Az  =  4.98  mm.  The  time  step  associated 
with  these  grid  dimensions  is  8.06  picoseconds.  The  difference  equations  which  are 
utilized  to  update  the  electric  and  magnetic  fields  in  a  leap  frog  manner  are  obtained 
from  Maxwell’s  equations  in  differential  form  by  approximating  spatial  and  temporal 
derivatives  with  center  differences.  Electric  and  magnetic  fields  are  discretized  on 
Yee  s  lattice,  which  places  them  such  that  they  are  interlaced  in  time  and  space  [1]. 
Initial  conditions,  which  are  imposed  at  the  beginning  of  the  simulation,  are  that  all 
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electric  and  magnetic  fields  start  at  zero.  At  each  subsequent  time  step,  fields  are 
imposed  at  appropriate  locations  in  order  to  simulate  the  excitation  source.  In  this 
problem,  the  source  is  implemented  at  the  wire  (connector)  connecting  the  bottom 
of  the  module  to  the  backplane  by  adjusting  the  electric  field  value  which  is  parallel 
with  and  coincides  with  the  wire  location  [25].  Assuming  that  the  wire  is  oriented 
in  the  2:  direction,  then  the  difference  equation  utilized  to  update  the  electric  field  at 
the  source  location  is 


k  +  1/2)  = 

[£  + 1/2,  /.  k  + 1/2)  -  - 1/2.  /.  k  + 1/2)) 

+  1/2.  >=  +  1/2)  -  -  1/2,  k  +  1/2))]  (4.1) 


where  n  is  the  temporal  index,  i,  j,  and  k  are  the  spatial  indices  in  the  x,  y,  and 
z  directions.  E^,  and  Hy  simply  represent  the  electric  and  magnetic  fields  in 
the  directions  corresponding  to  their  subscripts.  E^  represents  the  excitation  source 
and  is  chosen  to  be  a  Gaussian  pulse  modulated  by  a  sine  wave  in  time  and  hence 
contains  energy  in  a  band  of  frequencies  allowing  multiple  frequencies  to  be  analyzed 
simultaneously.  A  time  profile  for  Eg  would  have  the  form 


(4.2) 


with  typical  values  of  a  =  0.6  ns  and  /o  =  800  MHz.  cq  is  the  free  space  permittivity 
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and  a  is  the  conductivity  which  is  chosen  as 

_  Az 
ZoAxAy 

where  Zq  is  the  characteristic  impedance  of  the  line,  and  was  chosen  to  be  50f2  in 
order  to  match  the  coaxial  line  used  in  the  measurements.  This  equation  models 
a  matched  source  coaxial  cable  transmission  line  excitation.  Alternatively,  setting 
the  electric  field  to  have  a  predetermined  time  profile  is  equivalent  to  a  zero  source 
impedance  coaxial  transmission  line.  So  the  source  used  to  drive  this  configuration  is 
a  50Q  coaxial  cable  transmission  line  with  its  center  conductor  connected  to  the  mod¬ 
ule  and  the  outer  conductor  connected  to  the  backplane.  Boundary  conditions  are 
imposed  at  perfectly  conducting  surfaces  and  at  the  outer  boundary.  At  the  perfectly 
conducting  surfaces,  the  tangential  electric  field  is  set  to  zero.  Absorbing  boundary 
conditions  must  be  applied  at  the  outer  boundary  of  the  computational  domain  in 
order  to  simulate  unbounded  space  [67]- [75].  Second-order  Mur  [67]  absorbing  bound¬ 
ary  conditions  (ABC)  are  used  at  the  outer  boundary,  except  at  the  corners  and  edges 
where  first-order  ABCs  are  used.  Finally,  by  applying  the  difference  equations  ob¬ 
tained  from  Maxwell’s  equations  in  addition  to  the  appropriate  initial  and  boundary 
conditions,  the  electric  and  magnetic  fields  can  be  calculated  everywhere  within  the 
computational  domain  for  each  time  step.  The  FD-TD  simulation  is  executed  until 
steady  state  is  reached,  which  for  excitations  with  pulse  time  profiles  occurs  when  the 
fields  in  the  computational  domain  have  become  negligibly  small.  The  simulations 
executed  for  this  paper  required  from  2000  to  6000  time  steps  to  achieve  steady  state. 

The  calculation  of  total  radiated  power  and  radiation  impedance  is  done  using 
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time  domain  results  from  the  FD-TD  simulation.  First,  the  electric  and  magnetic 
fields  are  calculated  everywhere  in  the  computational  domain.  Both  the  voltage  and 
current  on  the  connector  are  recorded,  and  then  their  frequency  domain  amplitudes 
are  calculated,  whose  ratio  determines  the  radiation  impedance  [25].  The  current  on 
the  connector  Icon  is  calculated  using  Ampere’s  Law. 

hon{t)  =  H{t)  ■  di  (4.4) 

connector 

The  voltage  V^on  is  defined  as 

^backplane _ 

Vcon{t)=  E{t)-dl  (4.5) 

^module 

In  the  FD-TD  scheme  with  the  excitation  being  applied  to  the  electric  field  J5'z(z,  j,  k), 
the  time  domain  current  and  voltage  are  calculated  using 

/n+l/2  ^  [i/;+V2(i  +  1/2,  j,  k  -b  1/2)  -  _  1/2,  k  +  1/2)] 

-Aa;  +  1/2,  k  +  1/2)  -  -  1/2,  k  +  1/2)]  (4.6) 

=  -AzE^{i,j,  k  +  1/2)  (4.7) 

Finally,  Fourier  transforming  the  voltage  and  current  and  taking  their  ratio,  the 
radiation  impedance  Zr^d  is  obtained. 


^rad(/)=i?rad(/)+7Xad(/) 


Vconif) 

Iconif) 


(4.8) 
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This  radiation  impedance  is  at  the  connector  between  the  module  and  the  backplane. 

The  total  radiated  power  can  also  be  obtained  from  the  complex  amplitudes  of 
the  current  and  voltage. 


Prad(/)=2  Re(^co„(/)Kon(/)) 


(4.9) 


4.3.2  Prony’s  Method 


In  applying  the  FD-TD  technique  to  certain  of  the  more  resonant  modules-on-back- 
plane  configurations,  a  large  number  of  time  steps  (on  the  order  of  thousands)  are 
required  in  order  to  achieve  steady  state.  An  extrapolation  scheme  based  on  Prony’s 
method  is  used  in  conjunction  with  a  relatively  short  FD-TD  time  response  to  effi¬ 
ciently  obtain  the  late  time  response  [114].  This  extrapolation  method  was  applied 
to  the  time  response  of  the  current  and  voltage  on  the  connector.  The  FD-TD  simu¬ 
lation  is  executed  long  enough  such  that  the  remaining  time  response  approximately 
contains  only  a  few  decaying  complex  exponentials.  Prony’s  method  approximately 
determines  the  complex  amplitudes  and  complex  frequencies  of  these  exponentials. 
Once  the  complex  amplitudes  and  frequencies  are  determined  then  the  time  response 
is  extrapolated  to  the  desired  extent.  In  analyzing  the  more  resonant  modules-on- 
backplane  configurations,  the  FD-TD  simulation  would  typically  be  executed  for  3000 
time  steps  and  then  Prony’s  method  would  be  applied  to  extrapolate  the  next  3000 
time  steps. 
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The  time  response  is  assumed  to  have  the  following  form: 


M 


M 


fin]  =  finAt)  =  Yi  =  Y 


(4.10) 


m=l 


m=l 


where  n  is  the  number  of  the  time  step,  At  is  the  size  of  a  single  time  step,  M  is  the 
number  of  complex  exponentials  being  used  to  estimate  the  spectral  content  of  the 
time  response,  and  A^  and  are  the  complex  amplitudes  and  complex  frequencies 
of  the  exponentials.  The  time  response  f[n]  can  be  written  as  a  polynomial  where 


M 


f[n]  = 

m—l 


(4.11) 


and 


^m.  —  ^ 


(4.12) 


A  characteristic  polynomial  is  used,  which  is  defined  as 


M 


p(z)  =  n",,(z-2i..)  =  -Epi^ 

1=0 


M-l 


(4.13) 


where  po  =  —  1- 

By  rearranging  the  two  previous  equations,  the  following  expressions  can  be 
obtained. 


M  M  (  M  ^ 

Y  {PM}  =  Y 

m=l  m=l  V  0  > 


(4.14) 
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MM  M 

0  =  -  =  J2f[n-  l]pi  (4.15) 

/=0  m=l  l—O 


M 

/W  =  E/[«-^]P<  (4.16) 

i=\ 


The  coefficients  pi  to  the  characteristic  equation,  can  be  found  using  a  least 
squares  fit  to  the  time  response  f[n].  The  time  step  used  in  conjunction  with  Prony’s 
method  is  typically  an  integer  multiple  of  the  time  step  used  in  the  FD-TD  simulation. 
The  time  step  only  has  to  be  small  enough  to  capture  the  variation  of  the  time 
response,  or  more  exactly. 


-  TT  <  Im(s,„)At  <  TT  (4.17) 

Once  the  coefficients  are  determined,  the  roots  Zm  of  the  characteristic  equation 
can  be  calculated.  Another  constraint  on  the  time  step  is  that  At  should  not  be  too 
small,  since  that  will  result  in  all  the  roots  Zm  being  very  close  to  1.  When  a  number 
of  the  roots  are  too  close  to  1  then  the  solution  for  the  roots  becomes  numerically 
unreliable.  Hence  trying  to  use  too  many  complex  exponentials  or  too  small  a  time 
step  will  cause  numerical  problems. 

Finally,  the  complex  amplitudes  Am  can  be  found  using  the  following  linear 


system  of  equations. 
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where  n  is  chosen  to  be  the  last  M  evenly  spaced  samples  of  the  time  response. 


4.4  Measurement  Procedure 

Measurements  were  performed  using  an  HP  8510  Network  Analyzer.  The  network 
analyzer  was  connected  to  the  modules-on-backplane  configuration  to  be  tested  by  a 
test  port  cable  and  semi-rigid  coaxial  cable  [see  Figure  4.3].  While  the  measurements 
were  not  conducted  within  an  anechoic  chamber,  absorbers  were  utilized  to  reduce 
the  effects  of  the  surrounding  environment.  In  addition  to  the  room  effects,  the 
presence  of  the  test  port  cable  and  semi-rigid  coaxial  cable  on  the  backside  of  the 
backplane  could  have  caused  differences  with  the  FD-TD  predictions,  since  all  the 
cables  were  absent  in  the  FD-TD  simulations.  However,  since  the  backside  of  the 
backplane  is  shielded  by  the  backplane,  presumably  the  presence  of  the  cable  would 
not  significantly  affect  the  measurements.  In  each  of  the  measurements,  the  center 
conductor  of  the  coaxial  cable  was  soldered  to  the  module  under  test  such  that  the 
coaxial  cable  was  in  contact  with  at  least  one  of  the  reference  planes  in  the  module. 
5-parameter  data  (5ii)  was  recorded  using  this  set-up. 

The  radiation  impedance  is  deduced  from  these  measurements  by  performing 
the  following  steps.  The  reflection  coefficient,  F,  at  the  connector  between  the  module 
and  backplane  is  related  to  the  Su  data  which  is  obtained  at  the  end  of  the  test  port 
cable  by  the  following  phase  term. 
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Figure  4.3:  Measurement  setup. 

where  k  =  27r//ci,  /  is  the  frequency,  and  ci  is  the  speed  of  light  within  the  semi-rigid 
coaxial  cable  of  length  /.  Then  the  load  (radiation)  impedance  can  be  found  using 
this  relation. 


i  +  r 

(4.20) 

where  Zq  is  the  characteristic  impedance  of  the  coaxial  line.  The  characteristic 
impedance  of  the  line  used  in  our  study  was  50fi. 

The  following  two  configurations  were  examined:  monopoles  of  lengths  7.5  cm 
and  5  cm,  both  over  a  ground  plane  of  dimensions  30.48  cm  by  31.75  cm.  The  real 
and  imaginary  portions  of  the  radiation  impedance  of  the  7.5  cm  and  5  cm  monopoles 
over  a  finite  ground  plane  are  plotted  in  Figures  4.4  and  4.5,  respectively.  The  lowest 
resonances  of  these  structures,  which  correspond  to  purely  real  radiation  impedances, 
occur  at  0.913  GHz  and  1.393  GHz  for  the  7.5  cm  and  5  cm  monopoles,  respectively. 
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Figure  4.4:  Radiation  impedance  of  a  7.5  cm  monopole  over  a  ground  plane  of  di¬ 
mensions  30.48  cm  by  31.75  cm  as  a  function  of  frequency. 

The  radiation  impedances  at  resonance  are  34. 2Q  and  36.3f2  for  the  7.5  cm  and  5 
cm  monopoles,  respectively.  Theoretically,  the  resonant  frequencies  of  7.5  cm  and 
5  cm  monopoles  over  an  infinite  ground  plane  should  occur  at  1  GHz  and  1.5  GHz, 
respectively.  However,  due  to  finite  size  of  the  ground  planes,  the  resonant  frequencies 
are  slightly  shifted.  The  radiation  impedance  at  resonance  of  the  monopoles  over 
infinite  ground  planes  should  theoretically  be  about  36.5f2,  which  corresponds  well 


with  the  measured  data. 
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Figure  4.5:  Radiation  impedance  of  a  5  cm  monopole  over  a  ground  plane  of  dimen¬ 
sions  30.48  cm  by  31.75  cm  as  a  function  of  frequency. 
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4.5  Results  and  Discussion 

Measurements  are  obtained  for  single  module  configurations,  multiple  module  con¬ 
figurations  and  for  partially  enclosed  configurations.  FD-TD  numerical  results  have 
been  obtained  for  each  of  the  measured  configurations  and  fairly  good  agreement  has 
been  achieved. 

In  order  to  justify  the  simplified  model,  measurements  of  the  radiation  imped¬ 
ance  of  modules-on-backplane  configurations  with  plain  copper  modules  and  popu¬ 
lated  modules  are  compared.  The  populated  modules  were  standard  plug-in  printed 
circuit  boards  with  all  the  chips  and  devices  intact.  The  populated  module  was  ex¬ 
cited,  in  each  case,  by  the  center  conductor  of  a  coaxial  probe  which  was  connected 
to  a  reference  plane  in  the  printed  wire  board.  The  connection  was  accomplished 
by  soldering  the  center  conductor  to  a  pin  on  the  printed  wire  board  which  leads  to 
the  ground  plane  of  the  board.  The  outer  conductor  of  the  coaxial  line  was  soldered 
to  the  backplane  and  the  length  of  the  coaxial  probe  connector  extending  above  the 
backplane  in  this  case  and  in  all  subsequent  cases  is  1  cm.  Figures  4.6  and  4.7  show 
the  radiation  impedance  of  the  populated  and  plain  copper  modules  for  two  different 
feed  (source)  points.  The  feed  points  are  located  2.5  cm  away  from  the  center  of  the 
bottom  edge  of  the  module.  The  modules  have  dimensions  of  13.02  cm  in  width  and 
21.43  cm  in  height.  The  backplane  is  29.21  cm  by  30.48  cm.  The  module  is  centered 
on  the  backplane  and  the  width  of  the  module  is  parallel  to  the  29.21  cm  side  of 
the  backplane.  The  real  and  imaginary  portions  of  the  radiation  impedances  for  the 
populated  and  plain  copper  modules  show  good  correlation.  The  good  correlation 
between  the  radiation  impedances  for  the  populated  and  plain  copper  modules  at 
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Figure  4.6:  Radiation  impedance  of  a  populated  and  a  plain  copper  module  of  di¬ 
mension  13.01  cm  by  21.43  cm  over  a  finite  ground  plane  as  a  function  of  frequency 
with  the  source  connected  to  the  bottom  edge  of  the  left  side  of  the  module. 

both  feed  points  justifies  the  study  of  the  simplified  model. 

Figure  4.8  compares  the  radiation  impedance  as  a  function  of  frequency  obtained 
from  measurements  and  FD-TD  simulations  for  a  2.54  cm  wide  and  16.51  cm  high 
module  with  a  backplane  of  size  29.21  cm  by  30.48  cm.  In  this  case,  and  in  all 
subsequent  single  module  cases,  the  module  is  centered  on  the  backplane  with  the 
bottom  edge  of  the  module  parallel  to  the  29.21  cm  side  of  the  backplane.  The 
connector,  in  this  and  all  subsequent  cases,  is  attached  to  the  middle  of  the  bottom 
edge  of  the  module.  This  first  case  has  dimensions  atypical  of  a  standard  plug-in 
printed  circuit  board  module.  However,  this  case  was  chosen,  because  the  module 
is  thin  and  tall,  to  determine  if  the  measurements  and  FD-TD  predictions  would 
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Figure  4.7:  Radiation  impedance  of  a  populated  and  a  plain  copper  module  of  di¬ 
mension  13.01  cm  by  21.43  cm  over  a  finite  ground  plane  as  a  function  of  frequency 
with  the  source  connected  to  the  bottom  edge  of  the  right  side  of  the  module. 
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Figure  4.8:  Radiation  impedance  of  a  2.54  cm  by  16.51  cm  module  over  a  finite  ground 
plane  as  a  function  of  frequency. 

correspond  well  with  an  analytical  solution  of  a  monopole  over  an  infinite  ground 
plane.  The  first  resonance  occurs  at  375  MHz  (the  first  zero  crossing  in  the  reactance) 
which  is  somewhat  below  the  expected  quarter-wavelength  resonance  associated  with 
a  monopole  over  an  infinite  ground  plane  of  438  MHz.  This  discrepancy  is  primarily 
due  to  finite  nature  of  the  the  ground  plane.  The  next  resonance,  which  corresponds 
to  a  three  quarter-wavelength  resonance,  occurs  at  the  expected  frequency  of  1.125 
GHz,  which  is  three  times  higher  than  the  first  resonance.  Fairly  good  agreement 
was  achieved  for  the  radiation  impedance  between  the  FD-TD  prediction  and  the 
measurements. 

In  Figure  4.9,  Prony’s  method  is  demonstrated  for  the  previous  case  [Section 
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4.3.2].  The  time  response  of  the  current  on  the  connector  to  a  2.54  cm  by  16.51  cm 
module  over  a  29.21  cm  by  30.48  cm  ground  plane  is  used  as  a  sample  time  response. 
The  FD-TD  simulation  was  executed  for  3000  time  steps  and  Prony’s  method  was 
applied  using  only  the  first  1000  time  steps.  The  excitation  used  had  the  form  of 
Equation  4.2,  with  a  =  2  ns,  and  /o  =  1  GHz  with  a  peak  voltage  of  1  volt.  In 
applying  Prony’s  method,  nine  complex  exponentials  were  used  to  fit  the  data,  every 
tenth  time  step  was  utilized,  and  only  40  data  points  were  used.  From  Figure  4.9,  it 
is  apparent  that  the  Prony’s  method  extrapolation  fits  the  actual  FD-TD  simulation 
results  reasonably  well  but  shows  a  slight  phase  difference  as  time  progresses.  Note 
that  Figure  4.9  only  plots  the  time  response  after  the  point  at  which  Prony’s  method 
is  applied.  Prony’s  method  was  also  applied  to  the  voltage  on  the  connector  in 
order  to  obtain  the  radiation  impedance  using  the  Prony’s  method  extrapolated  time 
sequences.  Figure  4.10  plots  the  radiation  impedance  of  the  module,  and  illustrates: 
(1)  the  good  agreement  in  the  frequency  domain  between  the  long  FD-TD  simulation 
and  the  short  FD-TD  simulation  with  Prony’s  method  extrapolations,  and  (2)  the 
improvement  of  the  short  FD-TD  simulation  with  the  Prony’s  method  extrapolations. 

The  significance  of  searching  for  the  radiation  impedance  is  that  the  radiative 
properties  of  the  module  can  then  be  easily  incorporated  in  a  circuit  model  as  a 
lumped  element  with  a  complex  impedance  as  shown  in  Figure  4.11.  The  power 
radiated  for  the  configuration  can  be  expressed  in  terms  of  the  voltage  across  or  the 
current  through  the  radiator  model. 


Trad  = 


RrSid 


(4.21) 
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Figure  4.9;  Time  response  of  the  current  on  the  connector  to  a  2.54  cm  by  16.51  cm 
module  over  a  29.21  cm  by  30.49  cm  ground  plane  as  calculated  using  direct  FD-TD 
calculations  and  Prony’s  method  extrapolations. 
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Figure  4.10:  Comparison  of  radiation  impedance  for  a  2.54  cm  by  16.51  cm  module 
over  a  29.21  cm  by  30.49  cm  ground  plane  calculated  using  a  short  FD-TD  simula¬ 
tion,  a  FD-TD  simulation  with  Prony’s  method  extrapolations,  and  a  long  FD-TD 
simulation. 


4.5.  RESULTS  AND  DISCUSSION 


143 


CIRCUIT 

MODEL 


R 


RAD 


jx 

~T~ 


RADIATOR  MODEL 


Figure  4.11:  Characterization  of  radiative  properties  of  modules-on-backplane  as  cir¬ 
cuit  parameters. 


(4.22) 


So  when  a  voltage  source  is  used,  the  reactance  of  the  radiator  becomes  an  important 
factor  in  determining  the  amount  of  radiation  emitted. 

In  practice  the  radiated  power  is  of  greater  interest  than  the  radiation  imped¬ 
ance.  In  determining  the  radiated  power,  a  1  mV  voltage  source  is  utilized  (i.e.  the 
complex  amplitude  at  all  frequencies  was  1  mV),  with  zero  source  impedance.  Since 
the  power  distribution  to  the  reference  planes  is  designed  to  have  low  impedance,  a 
zero  source  impedance  voltage  source  is  a  reasonable  model  [113].  Figure  4.12  shows 
the  radiated  power  in  dBnW.  The  peak  radiation  occurs  at  the  first  resonance  around 
375  MHz  and  the  resulting  radiation  is  about  37  nW  for  a  1  mV  source. 

The  radiated  power  for  a  2.54  cm  wide  and  29.21  cm  high  module  with  a  back¬ 
plane  of  size  29.21  cm  by  30.48  cm  is  examined  in  Figure  4.13.  This  case,  which  also 
has  atypical  dimensions,  was  examined  to  determine  if  the  measurements  and  FD- 
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Figure  4.12:  Radiated  power  from  a  2.54  cm  by  16.51  cm  module  over  a  finite  ground 
plane  excited  by  a  1  mV  voltage  source  with  zero  impedance  as  a  function  of  frequency. 
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Figure  4.13:  Radiated  power  from  a  2.54  cm  by  29.21  cm  module  over  a  finite  ground 
plane  excited  by  a  1  mV  voltage  source  with  zero  impedance  as  a  function  of  frequency. 

TD  predictions  displayed  the  expected  trends  predicted  by  an  analytical  solution  of 
a  monopole  over  an  infinite  ground  plane.  Results  obtained  from  FD-TD  simulations 
and  from  measurements  are  plotted.  The  measured  data  agrees  fairly  well  with  the 
FD-TD  simulation  calculations.  The  first  resonance  occurs  at  240  MHz  which  again  is 
slightly  below  the  expected  quarter-wavelength  of  a  monopole  over  an  infinite  ground 
plane  of  248  MHz.  Again  resonances  are  observed  at  approximately  odd  multiples  of 
the  first  resonance,  which  is  consistent  with  a  monopole  over  an  infinite  ground  plane 
model.  As  expected,  more  peaks  in  the  radiation  resistance  and  zero  crossings  in  the 
radiation  reactance  are  apparent  compared  to  the  smaller  module. 


Figure  4.14  compares  the  radiated  power  obtained  from  measurements  and  from 
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Figure  4.14:  Radiated  power  from  a  13.01  cm  by  21.43  cm  module  over  a  finite  ground 
plane  by  a  1  mV  voltage  source  with  zero  impedance  as  a  function  of  frequency. 

FD-TD  simulations  for  a  13.01  cm  by  21.43  cm  module  with  a  backplane  of  size  29.21 
cm  by  30.48  cm.  This  module  represents  a  standard  plug-in  printed  wire  board 
module.  Due  to  the  more  squarish  dimensions,  the  module  will  no  longer  appear  to 
be  a  simple  monopole.  The  first  resonance  of  280  MHz  still  roughly  corresponds  to  the 
quarter-wavelength  monopole  over  a  ground  plane  resonance  of  334  MHz.  However, 
the  higher  resonances  are  more  difficult  to  explain,  as  resonances  across  the  module 
become  important  at  the  higher  frequencies.  The  measured  data  agree  fairly  well 
with  the  FD-TD  simulation  calculations.  Although  for  frequencies  above  0.7  GHz, 
there  is  a  3  to  5  dB  discrepancy  in  power  between  the  measurements  and  the  FD-TD 
predictions. 
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In  the  previous  cases,  single  modules  were  examined,  so  that  the  physics  of  the 
simpler  cases  could  be  understood.  In  examining  multiple  modules,  only  two  and 
three  module  configurations  were  considered.  The  modules  are  aligned  to  be  parallel 
to  each  other,  as  shown  in  Figure  4.2.  In  the  multiple  module  cases,  the  collection 
of  modules  are  centered  on  the  backplane,  with  the  bottom  edges  of  the  modules 
parallel  to  the  shorter  side  of  the  backplane.  The  modules  are  13.01  cm  by  21.43  cm 
in  size.  The  separation  between  the  modules  is  2.54  cm.  Only  one  of  the  multiple 
modules  is  excited  and  in  the  three  module  case  the  center  module  is  excited.  The 
unexcited  modules  are  grounded  at  the  center  of  their  bottom  edges  to  the  backplane 
by  a  1  cm  wire.  Figures  4.15  and  4.16  compare  the  radiated  powers  obtained  from 
measurements  and  from  FD-TD  simulations  for  two  and  three  13.01  cm  by  21.43  cm 
modules  with  a  backplane  of  size  29.21  cm  by  30.48  cm.  The  addition  of  a  second 
and  a  third  module  each  significantly  enhance  the  radiated  power  at  the  resonant 
frequencies  of  the  individual  configuration.  In  particular,  with  two  modules  there  is 
7  dB  enhancement,  and  with  three  modules  there  is  about  18  dB  enhancement  of 
radiation.  The  radiated  power  has  sharper  peaks  with  the  addition  of  the  second  and 
third  modules.  These  trends  are  consistent  with  the  fact  that  the  structure  resembles 
a  resonator  and  has  a  higher  Q  as  modules  are  added.  The  agreement  in  these  cases 
is  not  as  good  as  in  previous  cases,  although  the  general  shape  and  location  of  the 
resonances  of  the  measurements  and  FD-TD  predictions  are  the  same.  The  FD-TD 
predictions  tend  to  overestimate  the  radiation.  Possible  explanations  for  the  difference 
are:  First,  in  the  FD-TD  model  there  are  no  resistive  losses.  Second,  the  modules  in 
the  measurement  cases  may  not  have  been  precisely  aligned  in  that  the  modules  may 
not  have  been  perfectly  parallel  with  each  other,  which  would  tend  to  lower  the  Q 
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Figure  4.15:  Radiated  power  from  two  13.01  cm  by  21.43  cm  modules  over  a  finite 
ground  plane  excited  by  a  single  1  mV  voltage  source  with  zero  impedance  as  a 
function  of  frequency. 

factor. 

In  Figure  4.17,  the  effect  of  partially  enclosing  a  single  13.01  cm  by  21.43  cm 
module  is  examined.  The  structure  is  closed  on  the  bottom  and  four  sides  and  is  open 
on  top,  and  the  specific  dimensions  are  shown  in  Figure  4.18.  The  module  is  again 
excited  by  a  coaxial  probe  connected  to  the  center  of  the  bottom  edge  of  the  module. 
With  the  addition  of  the  enclosure,  the  radiated  power  of  this  configuration  has  sharp 
peaks  as  a  function  of  frequency,  with  a  peak  value  of  about  32  dBnW.  The  radiated 
power  predicted  by  the  FD-TD  simulation  shows  the  same  general  characteristics  as 
the  measured  radiated  power.  Similar  arguments  explaining  the  discrepancies  in  the 
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Figure  4.16.  Radiated  power  from  three  13.01  cm  by  21.43  cm  modules  over  a  finite 
ground  plane  excited  by  a  single  1  mV  voltage  source  with  zero  impedance  as  a 
function  of  frequency. 
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Figure  4.17:  Radiated  power  from  a  partially  enclosed  13.01  cm  by  21.43  cm  module 
over  a  finite  ground  plane  excited  by  a  1  mV  voltage  source  with  zero  impedance  as 
a  function  of  frequency. 

multiple  module  cases  can  be  applied  to  the  enclosed  configuration. 


4.5.  RESULTS  AND  DISCUSSION 


9.91  cm 


Figure  4.18:  Partially  enclosed  module-on-backplane  configuration. 
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4.6  Conclusions 

The  radiation  properties  of  modules-on-backplane  configurations  have  been  analyzed 
by  using  the  FD-TD  technique  and  by  taking  measurements.  The  FD-TD  technique 
was  chosen  for  its  flexibility  in  analyzing  this  class  of  potentially  complex  problems, 
for  which  neither  approximate  nor  analytical  solutions  are  applicable.  Fairly  good 
agreement  was  obtained  between  the  measurements  and  the  FD-TD  analysis.  The 
FD-TD  analysis  provides  radiation  impedances,  which  when  used  as  circuit  model 
parameters  can  be  utilized  to  estimate  overall  radiated  power.  The  FD-TD  analysis 
also  provides  the  capability  of  determining  near-fields  and  far-field  radiation  patterns. 

The  analysis  of  the  modules-on-backplane  configuration,  a  common  structure 
in  the  computer  environment,  has  shown  that  radiation  enhancement  occurs  in  the 
hundreds  of  megahertz.  The  addition  of  modules  or  an  enclosure  have  been  shown  to 
enhance  the  radiation  from  the  structure.  The  potential  for  increased  radiation  due  to 
the  modules-on-backplane  configuration  has  been  demonstrated  with  measurements 


and  numerical  simulations. 


Chapter  5 

Coupling  of  EM  Energy  by  a  Wire 
Penetrating  a  Metallic  Screen 

5.1  Introduction 


The  primary  motivation  for  studying  various  structures  which  compose  a  computer 
is  that  computers  are  now  operating  at  faster  speeds  and  consuming  more  power 
resulting  in  significantly  increased  levels  of  power  at  higher  frequencies.  Metallic 
enclosures  are  commonly  used  to  house  computer  systems  and  always  have  cables 
and  wires  exiting  from  them,  A  simplified  configuration  of  the  problem  of  interest  is 
shown  in  Figure  5.1.  Since  the  metallic  enclosure  functions  both  as  a  electromagnetic 
shield  and  as  a  resonator  structure,  the  most  important  sources  of  radiation  come 
from  any  holes  in  the  box  which  are  required  for  ventilation  and  from  any  cables 
exiting  from  the  box  which  are  required  for  power  and  communications.  In  particular 
for  cables,  the  radiation  sources  inside  the  enclosure  can  couple  to  a  wire  exiting  the 
enclosure,  and  the  wire  can  function  as  an  antenna  in  the  hundreds  of  MHz  frequency 
range.  For  these  reasons,  the  coupling  of  electromagnetic  energy  out  of  the  metallic 
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MODULES-ON-BACKPLANE  AS  SOURCE 


Figure  5.1:  A  simple  wire  exiting  a  metal  enclosure  configuration. 

enclosure  by  cables  is  of  great  interest.  The  objective  of  this  research  is  to  be  able  to 
quantify  the  amount  of  radiation  emitted  by  these  cables  as  a  function  of  frequency. 

In  analyzing  the  problem,  the  following  assumptions  and  simplifications  were 
used.  The  radiation  admittance  of  the  external  portion  of  the  wire  is  assumed  to 
be  simple  enough  such  that  it  can  be  adequately  modeled  by  a  set  of  of  series  RLC 
circuits  which  are  in  parallel  (lumped  element  model).  Essentially,  one  RLC  circuit 
is  used  to  match  each  resonance  present  in  the  wire’s  radiation  admittance.  An 
additional  assumption  is  that  the  radiation  mechanism  is  primarily  due  to  the  current 
on  the  wire  and  not  the  fields  traveling  through  the  aperture,  since  the  aperture  size 
is  only  being  modeled  in  a  simple  manner  and  will  depend  on  the  FD-TD  cell  size. 
Since  the  aperture  has  dimensions  much  smaller  than  those  of  the  external  wire,  then 
the  radiation  through  the  aperture  at  frequencies  where  the  external  wire  resonates 
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will  be  in  the  Rayleigh  region.  Hence,  the  radiation  through  the  aperture  should  be 
very  small.  The  ground  planes  and  wires  are  all  assumed  to  be  perfectly  conducting. 

The  motivation  for  using  the  lumped  element  model  within  the  FD-TD  tech¬ 
nique  instead  of  a  full  FD-TD  simulation  is  to  obtain  significant  computational  sav¬ 
ings.  In  particular,  the  lumped  element  model  replaces  the  model  of  the  region  outside 
the  enclosure,  thus  eliminating  the  need  to  include  the  outer  region  in  the  FD-TD 
simulation.  If  the  external  wire  is  significantly  large  or  larger  than  the  enclosure, 
then  the  reduction  in  the  computational  domain  size  can  be  significant.  In  analyzing 
a  completely  enclosed  metal  box,  the  FD-TD  technique  does  not  need  to  impose  any 
absorbing  boundary  conditions  since  the  metal  enclosure  naturally  imposes  conditions 
at  the  outer  boundary.  As  a  consequence  of  this,  the  FD-TD  technique  also  does  not 
need  the  typical  free  space  buffer  between  the  object  and  the  absorbing  boundary. 
However,  in  performing  a  simulation  of  an  enclosure  with  a  wire  exiting  it,  a  large 
free  space  region  surrounding  the  enclosure  and  the  wire  must  be  included  in  order  to 
facilitate  the  implementation  of  the  absorbing  boundary  condition.  For  example,  the 
computational  domain  required  to  model  a  one  meter  length  external  wire  attached  to 
a  0.75m  by  0.75m  by  0.15m  metal  enclosure  would  be,  on  the  order  of,  10  times  larger 
than  the  one  required  for  just  the  metal  enclosure.  The  computational  requirements 
of  the  lumped  element  model  are  a  mere  fraction  of  the  FD-TD  requirements,  on  the 
order  of  tens  of  storage  variables  compared  to  tens  of  thousands  of  storage  variables, 
and  similar  ratios  for  arithmetic  operations. 

In  applying  the  lumped  element  model,  the  original  wire  exiting  an  enclosure 
configuration  has  been  split,  requiring  both  the  internal  and  external  regions  to  be 
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solved.  The  external  region  can  be  solved  using  any  method  which  is  most  efficient 
or  convenient.  In  addition,  if  the  external  region  is  simple  enough,  the  radiation 
admittance  of  the  wire  can  be  approximated  using  analytical  solutions  for  wires  over 
infinite  ground  planes.  Hence,  the  external  problem  can  be  solved  first,  and  then 
the  internal  region  can  be  analyzed  using  the  FD-TD  technique  with  the  appropriate 
lumped  element  model.  Also  by  splitting  the  problem  into  two  regions,  the  probable 
problem  of  a  significant  difference  in  scale  between  features  in  the  two  regions  is 
naturally  solved,  as  each  region  can  be  solved  separately  using  whatever  scale  is  most 
appropriate  for  each  region.  This  natural  solution  obviates  the  need  for  using  variable 
size  grid  cells  within  one  computational  domain  in  order  to  accommodate  the  scale 
difference  in  the  external  and  internal  regions. 

In  Section  5.2,  the  models  used  to  examine  the  wire  penetrating  a  metal  enclo¬ 
sure  are  outlined.  Section  5.3  discusses  the  finite  difference-time  domain  technique 
briefly,  analytical  solutions  for  wires  over  a  ground  plane,  and  the  lumped  element 
model  used  in  conjunction  with  the  FD-TD  technique.  The  results  include  the  radia¬ 
tion  admittance  and  total  radiated  power  as  functions  of  frequency,  and  the  lumped 
element  FD-TD  simulation  results  are  compared  with  the  full  FD-TD  simulation  re¬ 
sults  in  Section  5.4.  Section  5.5  draws  overall  observations  regarding  the  accuracy  and 
usefulness  of  the  lumped  element  model  and  the  coupling  of  electromagnetic  energy 
by  wires  penetrating  a  metal  enclosure. 
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5.2  Problem  Configuration  and  Model 

Two  simplified  configurations  were  considered  in  analyzing  the  real  configuration. 
The  simplest  model  involves  a  wire  over  an  infinite  ground  plane  being  excited  by 
a  coaxial  line.  The  full  FD-TD  model  and  the  FD-TD  lumped  element  model  are 
shown  in  Figure  5.2.  This  first  model  is  used  to  validate  the  lumped  element  model, 
since  analytical  solutions  exist  for  a  complex  load  terminating  a  coaxial  transmission 
line.  This  model  will  also  provide  a  measure  for  how  well  the  lumped  element  can 
model  the  true  radiation  admittance.  The  next  model  involves  a  wire  connected  to 
an  infinite  metal  screen  exiting  through  a  small  aperture,  with  a  dipole  source  on  the 
“inner”  side  of  the  metal  screen.  Both  full  and  lumped  element  FD-TD  models  are 
shown  in  Figure  5.3  for  this  configuration. 
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Figure  5.2:  Full  FD-TD  model  and  lumped  element  FD-TD  model  for  a  wire  over  a 
ground  plane  fed  by  a  coaxial  line. 


5.2.  PROBLEM  CONFIGURATION  AND  MODEL 


159 


I 

DIPOLE 


INFINITE 

CONDUCTING 

SCREEN 

SMALL  APERTURE 


WIRE 


INFINITE 

CONDUCTING 

SCREEN 


t 

DIPOLE 


LUMPED 

ELEMENT 


Figure  5.3:  Full  FD-TD  model  and  lumped  element  FD-TD  model  for  a  wire  exiting 
a  small  aperture  in  a  metal  screen  excited  by  a  dipole. 
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5.3  Method  of  Solution 

5.3.1  Finite  Difference  -  Time  Domain  Technique 

The  finite  difference-time  domain  (FD-TD)  technique  was  used  to  analyze  this  prob¬ 
lem  [1]-[12].  This  method  was  chosen  because  it  has  the  capability  of  analyzing 
multiple  frequencies  with  a  single  simulation,  it  is  well  suited  for  low  loss  models  and 
rectilinear  geometries,  and  it  is  a  proven  method,  as  it  has  been  previously  used  to 
successfully  analyze  radiation  from  a  heatsink  configuration  [31]-[33]  and  radiation 
from  modules-on-backplane  configurations  [34], [35]. 

A  rectangular  three-dimensional  grid  was  used  as  the  computational  domain. 
Spatially,  the  computational  domain  is  discretized  such  that  the  largest  grid  dimen¬ 
sion  is  less  than  about  a  twentieth  of  the  shortest  wavelength  of  interest.  The  time 
step  must  be  approximately  less  than  the  smallest  grid  dimension  divided  by  the 
speed  of  light  in  order  to  guarantee  the  stability  of  the  algorithm.  The  difference 
equations  which  are  utilized  to  update  the  electric  and  magnetic  fields  are  obtained 
from  Maxwell’s  equations  in  differential  form  by  approximating  spatial  and  temporal 
derivatives  with  center  differences.  Electric  and  magnetic  fields  are  discretized  on 
Yee’s  lattice  [1].  Initial  conditions,  which  are  imposed  at  the  beginning  of  the  simula¬ 
tion,  are  that  all  electric  and  magnetic  fields  start  at  zero.  At  each  subsequent  time 
step,  fields  are  imposed  at  appropriate  locations  in  order  to  simulate  the  excitation 
source.  The  imposed  electric  field  has  a  time  dependence  of  a  Gaussian  pulse  modu¬ 
lated  by  a  sine  wave,  which  contains  energy  in  a  band  of  frequencies  allowing  multiple 
frequencies  to  be  analyzed  simultaneously.  Boundary  conditions  are  imposed  at  per- 
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fectly  conducting  surfaces  and  at  the  outer  boundary.  At  the  perfectly  conducting 
surfaces,  the  tangential  electric  field  is  set  to  zero.  Absorbing  boundary  conditions 
must  be  applied  at  the  outer  boundary  of  the  computational  domain  in  order  to  simu¬ 
late  unbounded  space.  Second-order  absorbing  boundary  conditions  (ABC)  are  used 
at  the  outer  boundary,  except  at  the  corners  and  edges  where  first-order  ABCs  are 
used  [67]-[75].  Finally,  by  applying  the  difference  equations  obtained  from  Maxwell’s 
equations  in  addition  to  the  appropriate  initial  and  boundary  conditions,  the  electric 
and  magnetic  fields  can  be  calculated  everywhere  within  the  computational  domain 
for  each  time  step. 

The  calculation  of  total  power  radiated  is  done  using  time  domain  results  from 
the  FD-TD  simulation.  First,  the  electric  and  magnetic  fields  are  calculated  every¬ 
where  in  the  computational  domain  and  then  the  current  and  voltage  on  the  wire  are 
recorded  at  the  aperture  of  the  screen.  The  current  on  the  probe,  which  models  the 
connector,  is  calculated  using  Ampere’s  Law. 

7aperture(0  =  ^  (^-1) 

probe 

and  the  voltage  on  the  probe  is 

router  conductor _ 

Aperture  (^)  =  /  E{f)  '  (R  (5-2) 

«/ center  conductor 


Fourier  transforming  the  voltage  and  current  to  obtain  their  complex  amplitudes. 
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their  product  give  the  power  radiated. 

Piadif)  —  2  1^6  (  ■^aperture(/)^perture(/)^  (^•^) 

The  radiation  admittance  can  also  be  calculated  from  the  frequency  domain  voltage 
and  current  using  [25] 


YUf)  = 


-^aperture  (/) 
Aperture  (/) 


(5.4) 


5.3.2  Analytical  Solution  for  Wire  over  a  Ground  Plane 


In  applying  the  lumped  element  model  discussed  in  the  next  section,  it  was  necessary 
to  obtain  the  radiation  admittance  for  the  wire  exiting  from  the  enclosure.  The 
simplest  wire  model  involves  a  wire  which  is  normal  to  an  infinite  ground  plane.  This 
configuration  has  analytic  solutions  for  the  radiation  admittance. 

The  analytical  expression  for  the  radiation  impedance  of  wire  over  an  infinite 
ground  plane  is  as  follows  [116]-[118]: 


Y  — 

J.  yfiYe  - 


^wire  iXy, 


(6.5) 


~  +0.5sin(*:Q(Si(2*:<)  -  2Si():Q] 


+0.5cos(H)[7  +  ln(H/2)  +  Ci{2k£)  —  2Ci(fc^)]) 


(5.6) 
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^wire  =  — ^7777^(2Si(A:^)  +  cos(H)[2Si(A:£)-Si(2M)] 

OTTSin  [Kt/Z) 

-  siniM)  [2Ci{ke)  -  Ci(2M)  -  Ci(2A:aV£)])  (5.7) 

where  k  is  the  wavenumber,  £  is  the  imaged  length  (or  twice  the  true  length  1)  of  the 
wire,  a  is  the  radius  of  the  wire,  770  =  377f2  is  the  free  space  impedance,  7  =  0.5772 
and  Si  and  Ci  are  sine  and  cosine  integrals. 

Figure  5.4  compares  the  theoretical  radiation  admittance  of  a  wire  with  a  length 
of  Z  =  0.5  m  and  a  radius  of  a  =  2  mm  with  the  radiation  admittance  obtained  from 
a  full  FD-TD  simulation.  The  figure  shows  good  agreement,  and  demonstrates  the 
usefulness  of  the  theoretical  solution  for  the  radiation  admittance  of  a  wire. 


5.3.3  Lumped  Element  Model 

The  implementation  of  single  Rs,  Ls  and  Cs  has  been  done  previously  [22], [23],  how¬ 
ever  the  modeling  of  a  collection  of  these  elements  at  a  single  node  is  new.  The 
lumped  element  model  used  is  shown  in  Figure  5.5  and  consists  of  a  set  of  N  series 
RLC  circuits  in  parallel.  The  total  admittance  has  the  form 

y  =  ^=j2Ym  =  G  +  iS  '  (5.8) 

^  m=l 

where  G  and  S  are  the  conductance  and  susceptance  respectively,  and  the  admittance 
of  each  individual  RLC  circuit  is 


_ S/Lm _ 

+  sRrji/ Lm  +  I/{LffiCfn} 


i;^(s)  = 


(5.9) 


ADMITTANCE  (mhos) 
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Figure  5.4:  Comparison  of  theoretical  vs.  full  FD-TD  simulation  of  radiation  admit¬ 
tance  of  a  wire  over  an  infinite  ground  plane  with  a  length  of  0.5  m  and  a  radius  2 


mm. 
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Figure  5.5:  Lumped  element  model. 


where  s  —  — i27r/. 

The  circuit  elements  for  a  single  RLC  circuit  are  chosen  to  match  one  resonance 
in  the  wire’s  radiation  admittance. 


Rm  = 


Qmax 
^  m 


Rn 


2jrA/„ 


(5.10) 


(5.11) 


Crr.  = 


(27r/„)2L. 


(5.12) 


where  is  the  peak  of  the  conductance  at  resonance,  A/„i  is  the  half-power 

bandwidth  of  the  resonance  in  Hz,  and  fm  is  the  resonant  frequency  in  Hz. 
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The  transfer  function  relating  the  voltage  to  the  current  in  the  s  domain  is 


LjjiS  +  Rm  "h 


(5.13) 


where 


(5.14) 

m—1 

The  equivalent  equation  in  the  time  domain  is 

V{t)  =  +  Rm-^  [  Im{t)  (5.15) 

This  equation  needs  to  be  discretized,  and  backwards  differences  are  applied  to  insure 
stability.  The  resulting  equation  is 

(rn+1/2  _  rn-l/2\  \4.  n 

-=! - (5.16) 

The  update  equation  for  the  current  is 

Tn+l/2  _  ^  Y  ,  Lm/^t  jn-i/2  _  A  W-i/2  /r  jy'v 

LJAt  +  R^  Lm/At  +  Rm  + 

In  order  to  incorporate  the  above  difference  equation  which  is  in  terms  of  volt¬ 
ages  and  currents  into  the  FD-TD  scheme,  the  voltages  and  currents  need  to  be 
converted  to  electric  field  and  current  densities. 


jn+i/2  ^ 


(5.18) 
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^n+l/2  _ 


yn+l  _  yn  £n+l 

- -  =  ^ — - 2-Az 


(5.19) 


The  standard  difference  equation  for  updating  is 

,  .  E^*%j,k  +  l/2)-E^(i,j,k  +  l/2) 

- - 

+  1/2,  j,  A  +  1/2)  -  H;+'/^i  -l/2,j,k  +  1/2) 

Ax 

+  1/2,  k  +  1/2)  -  -  1/2,  k  +  1/2) 

Ay 

+ 1/2)  (5.20) 


Combining  the  expressions  for  the  voltage  in  terms  of  the  electric  field,  and  the 
current  in  terms  of  the  current  density,  the  following  equation  for  updating  E^  is 
obtained. 


B?+‘(i,j.i  +  l/2) 


Az 

.  At  2  Ax  Ay  (Lm/ At  +  Rm) 


E^{i,j,k  + 1/2) 


2AxAy{LjAt  +  R^) 


+  1/2,  j,  k  +  1/2)  -  //»+V^(z  -  1/2,  j,  k  +  1/2) 
Aa; 

+  1/2,  k  +  1/2)  -  j  -  1/2,  k  +  1/2) 

Ay 


N 


LJAt 


^ml H-  Rfri 


rn-1/2 

fim 


At/Cr, 


Lm/ +  Rui 


(5.21) 
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Resonance 

/o  (GHz) 

Gmax  (mhos) 

A/  (GHz) 

1 

0.1465 

2.9096  X  10-2 

0.0210 

2 

0.4420 

2.0300  X  10-2 

0.0355 

3 

0.7410 

1.7582  X  10-2 

0.0438 

4 

1.0405 

1.6129  X  10-2 

0.0497 

5 

1.3400 

1.5184  X  10-2 

0.0543 

6 

1.6395 

1.4500  X  10-2 

0.0580 

7 

1.9395 

1.3976  X  10-2 

0.0609 

Table  5.1:  Theoretical  resonances  for  a  wire  of  length  0.5  m  and  radius  of  2  mm  over 
an  infinite  ground  plane. 


5.4  Results  and  Discussion 


The  two  configurations  examined  include  a  wire  over  an  infinite  ground  plane  fed 
by  a  coaxial  line  [Figure  5.2]  and  a  grounded  wire  exiting  an  infinite  metal  screen 
through  a  small  aperture  [Figure  5.3].  For  both  cases,  an  identical  wire  configuration 
is  used  on  the  external  (right-hand)  side  of  each  model.  The  external  wire  extends  in 
a  direction  normal  to  the  infinite  ground  plane  and  has  a  length  of  /  =  0.5  m  and  a 
radius  of  a  =  2  mm. 

In  order  to  obtain  a  lumped  element  approximation  to  the  radiation  admit¬ 
tance  of  this  wire,  it  is  necessary  to  characterize  its  radiation  admittance.  Table 
5.1  describes  the  peaks  in  the  radiation  admittance  of  the  wire  antenna.  From  the 
theoretical  radiation  admittance,  the  following  RLC  lumped  element  model  was  con¬ 
structed.  Table  5.2  contains  the  values  of  the  circuit  elements  used  to  approximate 
the  wire  antenna’s  admittance. 


Figure  5.6  shows  the  analytical  solution  for  the  radiation  admittance  of  the  wire 
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Resonance 

R  {D) 

L(H) 

C(F) 

1 

34.36 

2.6047  X  lO-’^ 

4.5310  X  10-^=^ 

2 

49.26 

2.2084  X  10-'^ 

5.8708  X  10-13 

3 

56.87 

2.0667  X  10-^ 

2.2321  X  10-13 

4 

62.00 

1.9854  X  10-^ 

1.1784  X  10-13 

5 

65.85 

1.9303  X  10-^ 

7.3079  X  lO-i'i 

6 

68.96 

1.8924  X  10"^ 

4.9795  X  10-1^ 

7 

71.55 

1.8699  X  10-^ 

3.6011  X  10-1^ 

Table  5.2:  RLC  circuit  parameters  used  to  model  the  admittance  of  a  wire  of  length 
0.5  m  and  radius  2  mm  over  an  infinite  ground  plane. 

and  the  admittance  of  the  single  RLC  circuit.  The  match  appears  to  be  reasonably 
good  up  to  about  240  MHz.  By  using  2  RLC  circuits,  a  reasonable  match  is  obtained 
between  the  radiation  admittance  of  the  wire  and  the  RLC  circuit,  up  to  570  MHz 
as  shown  in  Figure  5.7.  Figure  5.8  shows  a  reasonable  match  up  to  860  MHz  between 
the  analytical  solution  for  the  radiation  admittance  of  the  wire  and  the  admittance 
of  the  three  parallel  RLC  circuits.  Finally,  Figure  5.9  shows  a  reasonable  match  up 
to  2000  MHz  between  the  analytical  solution  for  the  radiation  admittance  of  the  wire 
and  the  admittance  of  the  seven  parallel  RLC  circuits.  Figures  5.6  to  5.9  illustrate 
that  the  radiation  admittance  of  an  antenna  can  be  approximated  over  any  frequency 
range  with  the  use  of  additional  RLC  circuit  elements. 

The  simplest  configuration  which  is  shown  in  Figure  5.2  is  considered  first.  In 
this  case,  the  coaxial  line  is  excited  some  distance  below  the  ground  plane.  The  coaxial 
line  cross  section  is  shown  in  Figure  5.10.  A  cut  away  view  of  the  coaxial  line  is  shown 
in  Figure  5.11.  The  excitation  source  is  a  voltage  source,  14(0)  backed  by  a  perfect 
conductor  for  the  early  time  steps.  Once  the  wave  has  been  launched,  the  conductor 
is  replaced  by  a  first-order  absorbing  boundary  condition.  This  implementation  of 
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Figure  5.6:  Radiation  admittance  of  a  wire  of  length  0.5  m  and  radius  2  mm  over 
a  ground  plane  vs.  the  admittance  of  a  lumped  element  consisting  of  a  single  RLC 
circuit. 
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Figure  5.7:  Radiation  admittance  of  a  wire  of  length  0.5  m  and  radius  2  mm  over  a 
ground  plane  vs.  the  admittance  of  a  lumped  element  consisting  of  two  RLC  circuits. 
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Figure  5.8:  Radiation  admittance  of  a  wire  of  length  0.5  m  and  radius  2  mm  over 
a  ground  plane  vs.  the  admittance  of  a  lumped  element  consisting  of  three  RLC 
circuits. 
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Figure  5.9:  Radiation  admittance  of  a  wire  of  length  0.5  m  and  radius  2  mm  over 
a  ground  plane  vs.  the  admittance  of  a  lumped  element  consisting  of  seven  RLC 
circuits. 
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5  cm 


Figure  5.10;  Cross  section  of  the  coaxial  line. 

the  excitation  source  is  the  equivalent  of  launching  a  forward  traveling  wave  with 
V+(t)  =  Vo(t)  on  an  infinite  transmission  line,  or  a  matched  source  with  Vs(t)  =  2Vo{t). 

Figure  5.12  compares  the  lumped  element  FD-TD  results  and  the  circuit  solu¬ 
tions  for  the  input  admittance  at  the  termination  of  the  coaxial  line  for  seven  RLC 
lumped  elements.  The  figure  shows  near  perfect  agreement  between  the  FD-TD 
lumped  element  solution  and  the  circuit  solution  for  the  input  admittance  at  the  ter¬ 
mination  of  the  coaxial  line.  For  Figure  5.12,  the  cell  dimensions  used  in  the  FD-TD 
simulation  were  Ax  =  Ay  =  Az  =  0.0125m.  The  characteristic  impedance  of  the 
coaxial  line  was  Zq  =  141.5  f2.  An  additional  test  of  the  lumped  element  model  was 
performed  by  comparing  the  power  dissipated  in  the  lumped  element  calculated  by 
two  different  methods.  The  first  method  involved  simply  using  Equation  5.3,  which 
uses  the  current  and  voltage  at  the  lumped  element  to  determined  the  dissipated 
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Figure  5.11:  Cut  away  view  of  coaxial  line. 


power.  The  second  method  used  the  lumped  element  FD-TD  calculated  admittance, 


the  known  excitation  profile  on  the  transmission  line,  and  the  characteristic 


impedance  of  the  transmission  line  Zq,  in  conjunction  with  transmission  line  theory 


to  determined  the  dissipated  power.  Figure  5.13  shows  the  comparison  between  the 


two  methods  for  calculating  the  dissipated  power.  The  source  is  normalized  such 


that  the  complex  amplitude  of  the  exciting  matched  voltage  source.  Vs,  is  2  volts  at 


all  frequencies.  The  excellent  agreement  illustrates  that  the  lumped  element  model 


exhibits  the  correct  behavior  with  respect  to  the  fields  and  power. 


For  the  given  grid  cell  size,  the  coaxial  line  was  atypically  large.  This  large 


coaxial  line  was  examined  so  that  a  comparison  with  the  full  FD-TD  simulation,  which 


requires  the  larger  cell  size  in  order  to  have  a  tractable  size  computational  domain. 


would  be  feasible.  Figure  5.14  compares  the  radiation  admittance  obtained  using  a 
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Figure  5.12:  Input  admittance  at  the  termination  of  the  coaxial  line.  Lumped  element 
FD-TD  with  seven  RLC  vs.  circuit  solution. 
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Figure  5.13:  Radiated  power  of  a  coaxial  line  fed  wire  of  length  0.5  m  and  radius  2 
mm  for  the  lumped  element  FD-TD  with  seven  RLC  circuits  calculated  using  the 
admittance  vs.  calculated  using  the  current  and  voltage. 
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Figure  5.14:  Radiation  admittance  of  coaxial  line  fed  wire  of  length  0.5  m  and  radius  2 
mm  obtained  using  lumped  element  FD-TD  with  seven  RLC  circuits  vs.  full  FD-TD. 

lumped  element  FD-TD  simulation  and  a  full  FD-TD  simulation.  The  agreement 
here  appears  to  be  relatively  good,  as  all  the  peaks  in  the  radiation  conductance  are 
matched  in  magnitude  and  in  frequency. 

In  Figure  5.15,  the  radiated  power  is  compared  between  the  full  FD-TD  simula¬ 
tion  and  lumped  element  FD-TD  simulation.  In  both  cases,  the  source  is  normalized 
such  that  the  complex  amplitude  of  the  exciting  matched  voltage  source,  14,  is  2  volts 
at  all  frequencies.  For  reference,  the  power  dissipated  in  a  matched  load  would  be 
—24.5  dBW  for  this  transmission  line  and  source.  Figure  5.15  shows  that  a  significant 
portion  of  the  power  delivered  to  the  lumped  element  is  radiated  (only  1  to  2  dB  less 
then  the  maximum  possible)  at  frequencies  at  which  the  radiation  admittance  peaks 
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and  is  purely  real.  The  overall  levels  and  shape  correspond  fairly  well,  however,  the 
locations  of  the  peaks  show  a  slight  difference.  The  radiated  power  is  determined  by 
the  mismatch  between  the  coaxial  line  admittance  of  Yq  =  7.07  mi7~^,  and  the  input 
admittance  at  the  termination  of  the  coaxial  line.  In  the  lumped  element  FD-TD 
simulation,  the  terminating  admittance  is  the  lumped  element.  In  the  full  FD-TD 
simulation,  the  terminating  admittance  is  the  radiation  admittance  of  the  wire  over 
the  infinite  ground  plane.  The  mismatch  for  these  cases  is  minimized  when  the  ter¬ 
minating  admittance  is  closest  to  being  the  complex  conjugate  of  Yq.  So  for  these  two 
cases,  the  radiation  is  maximized  when  the  terminating  admittance  is  real  and  at  a 
maximum.  The  maxima  in  radiated  power  for  the  full  FD-TD  simulation  occur  at 
higher  frequencies  than  the  the  maxima  for  the  lumped  element  FD-TD  because  the 
zero  crossings  in  reactance  associated  with  the  peaks  in  the  conductance  for  the  full 
FD-TD  simulation  also  occur  at  slightly  higher  frequencies  than  those  associated  with 
the  lumped  element  FD-TD  simulation.  The  maxima  in  radiated  power  calculated 
by  the  two  methods  are  very  similar  in  magnitude  because  the  peaks  in  the  radiation 
admittance  are  also  very  similar  in  magnitude.  The  minima  in  the  radiated  power 
occur  at  the  zero  crossings  in  the  reactance  associated  with  the  minima  in  the  conduc¬ 
tance.  Since  at  these  frequencies,  the  lumped  element  FD-TD  calculated  conductance 
is  much  lower  than  the  full  FD-TD  calculated  conductance,  the  mismatch  with  the 
transmission  line  admittance  is  greater,  which  results  in  lower  radiated  power  for  the 
lumped  element  FD-TD  simulation.  Better  correlation  can  be  obtained  by  improving 
the  lumped  element  model  to  better  match  the  radiation  admittance  of  the  wire. 


The  next  configurations  involve  an  infinite  metal  screen,  and  their  geometry 


180 


CHAPTER  5.  EM  COUPLING  BY  WIRE  PENETRATING  SCREEN 


Figure  5.15:  Radiated  power  of  a  coaxial  line  fed  wire  obtained  using  lumped  element 
FD-TD  with  seven  RLC  circuits  vs.  full  FD-TD. 
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Figure  5.16:  Dimensions  of  wire  penetrating  a  metal  screen  configuration. 

is  shown  in  Figure  5.16.  For  configuration  A,  the  length  of  the  external  wire  is 
I  =  50  cm.  The  internal  wire  loop  has  height,  h  =  12.5  cm,  and  width  w  =  12.5 
cm.  The  exciting  dipole  is  d  =  25  cm  displaced  from  the  metal  screen.  The  cell 
dimensions  for  FD-TD  simulations  for  this  configuration  were  Ax  =  Ay  =  0.0125m 
and  Az  —  0.025m.  Figure  5.17  compares  the  radiation  admittance  obtained  using 
the  lumped  element  FD-TD  model  and  using  a  full  FD-TD  model.  The  agreement 
is  fairly  good  and  comparable  to  the  agreement  seen  earlier  between  the  theoretical 
radiation  admittance  of  the  wire  and  its  lumped  element  approximation. 

In  Figure  5.18,  the  radiated  power  is  compared  between  the  full  FD-TD  sim¬ 
ulation  and  lumped  element  FD-TD  simulation  for  configuration  A.  The  source  is 
normalized  such  that  the  power  radiated  from  the  exciting  dipole  in  free  space  is 
constant  for  all  frequencies.  Again  the  overall  levels  and  shape  correspond  fairly  well. 
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Figure  5.17:  Radiation  admittance  using  a  lumped  element  FD-TD  simulation  with 
three  RLC  circuits  and  using  a  full  FD-TD  simulation  of  the  infinite  metal  screen 
configuration  A. 
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However,  the  locations  of  the  peaks  show  a  difference.  The  maximum  coupling  occurs 
near  100  MHz  and  the  amount  of  power  coupled  is  only  2  to  3  dB  less  than  the  total 
power  radiated  by  the  dipole  in  free  space.  Note  that  the  presence  of  the  metal  screen 
and  wires  will  alter  the  total  power  radiated  by  the  dipole.  The  radiated  power  from 
the  external  wire  will  be  determined  by  two  factors.  The  first  factor  is  the  coupling 
between  the  exciting  dipole  and  the  internal  wire  loop.  The  second  factor  is  the  mis¬ 
match  between  the  radiation  impedance  of  the  internal  wire  loop,  and  the  radiation 
(lumped  element)  impedance  of  the  external  wire.  Figure  5.19  shows  the  radiation 
admittance  of  the  internal  wire  loop  for  configuration  A.  Figures  5.19  and  5.17  show 
that  the  frequencies  at  which  the  two  admittances  are  closest  to  being  conjugates  of 
each  other,  approximately  coincide  with  the  frequences  at  which  the  radiated  power 
peaks. 

For  configuration  B  (Figure  5.16),  the  length  of  the  external  wire  is  again  Z  =  50 
cm.  However,  the  internal  wire  loop  has  height,  =  10  cm,  and  width  w  =  7.5 
cm.  The  exciting  dipole  is  d  =  25  cm  displaced  from  the  metal  screen.  The  same 
cell  dimensions  were  used  for  configuration  B  as  in  configuration  A.  The  radiation 
admittance  obtained  using  the  lumped  element  FD-TD  model  and  using  a  full  FD- 
TD  model  are  nearly  identical  to  those  shown  in  Figure  5.17,  which  is  expected  as 
the  external  wire  configuration  was  unchanged.  In  Figure  5.20,  the  radiated  power  is 
compared  between  the  full  FD-TD  simulation  and  lumped  element  FD-TD  simulation 
for  configuration  B.  Again  the  source  is  normalized  such  that  the  power  radiated  from 
the  exciting  dipole  in  free  space  is  constant  for  all  frequencies.  The  maximum  coupling 
occurs  again  near  100  MHz  and  the  amount  of  power  coupled  is  5  to  6  dB  less  than 
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Figure  5.18:  Radiated  power  using  a  lumped  element  FD-TD  simulation  with  three 
RLC  circuits  and  using  a  full  FD-TD  simulation  of  the  infinite  metal  screen  configu¬ 
ration  A. 
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Figure  5.20:  Radiated  power  using  a  lumped  element  FD-TD  simulation  with  three 
RLC  circuits  and  using  a  full  FD-TD  simulation  of  the  infinite  metal  screen  configu¬ 
ration  B. 

the  total  power  radiated  by  the  dipole  in  free  space.  The  reason  for  the  lower  coupling 
in  configuration  B  compared  to  configuration  A  is  simply  that  a  better  admittance 
match  occurs  for  configuration  A.  The  overall  levels  and  shape  correspond  fairly  well. 
Figure  5.21  shows  the  radiation  admittance  of  the  internal  wire  loop.  Figures  5.21 
and  5.17  show  that  the  frequencies  at  which  the  two  admittances  are  closest  to  being 
conjugates  of  each  other,  approximately  coincide  with  the  frequencies  at  which  the 
radiated  power  peaks. 

The  differences  between  the  lumped  element  FD-TD  simulation  and  the  full  FD- 
TD  simulation  can  be  attributed  to  the  inaccuracies  of  the  lumped  element  model. 
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In  particular,  since  the  radiation  susceptance  has  greater  variation  with  increasing 
frequency,  the  error  in  the  radiated  power,  as  seen  in  Figures  5.18  and  5.20,  is  expected 
to  increase  with  frequency  as  well. 
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5.5  Conclusions 


Radiation  levels  associated  with  the  two  configurations  were  quantified.  Significant 
coupling  of  electromagnetic  energy  occurs  when  the  admittance  of  the  external  wire 
is  matched  by  the  admittance  associated  with  internal  side.  In  the  coaxial  line  con¬ 
figuration,  the  characteristic  admittance  of  the  coaxial  line  matches  the  peaks  in  the 
external  wire’s  radiation  admittance  well,  and  results  in  radiated  power  only  1  to  2  dB 
less  than  the  maximum  possible.  For  the  infinite  metal  screen  configurations,  fairly 
significant  levels  of  radiation  occur  at  frequencies  where  the  radiation  admittance  of 
the  inner  wire  loop  match  the  radiation  admittance  of  the  external  wire  antenna.  For 
configuration  A,  the  maximum  coupling  is  only  3dB  less  than  the  power  radiated  by 
the  dipole  in  free  space,  and  for  configuration  B,  the  maximum  coupling  is  only  6  dB 
less. 

A  “smart”  cell  in  the  FD-TD  scheme,  which  can  model  the  radiation  admit¬ 
tance  of  a  wire  antenna,  has  been  developed  and  tested.  Analytical  solutions  for  the 
radiation  admittance  of  a  straight  wire  over  a  ground  plane  have  been  recited  and 
approximated  using  the  lumped  element  model,  which  consists  of  a  set  of  series  RLC 
circuits  in  parallel.  The  admittance  of  the  lumped  element  implemented  in  the  FD- 
TD  simulations  match  the  circuit  solution  admittance  extremely  well.  Both  a  wire 
over  a  ground  plane  excited  by  a  coaxial  line,  and  a  wire  penetrating  a  metal  screen 
excited  by  a  dipole,  were  examined.  Fairly  good  correlation  between  lumped  element 
FD-TD  simulations  and  full  FD-TD  simulations  have  been  obtained  for  each  of  the 
configurations  for  both  radiation  admittances  and  radiated  powers. 
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One  clear  limitation  of  the  use  of  lumped  elements  in  modeling  radiation  admit¬ 
tances  is  that  the  results  will  only  be  as  accurate  as  the  lumped  element  model  of  the 
radiation  admittance.  The  lumped  element  model  considered  here,  most  accurately 
modeled  the  peaks  in  the  radiation  conductance,  and  hence  produced  the  most  accu¬ 
rate  results  at  those  frequencies.  Correspondingly,  the  lumped  element  FD-TD  model 
had  increasing  errors  as  the  lumped  element  admittance  varied  from  the  theoretical 
radiation  admittance.  In  the  future,  if  methods  for  designing  passive  RLC  circuits 
to  more  accurately  approximate  arbitrary  admittances  as  a  function  of  frequency  are 
developed,  then  the  lumped  element  model  could  be  correspondingly  more  accurate. 


Chapter  6 


EM  Radiation  from  a  Microstrip 
Patch  Antenna  on  a  Non-Uniform 
Substrate 


6.1  Introduction 


This  chapter  investigates  microstrip  patch  antennas  on  non-uniform  substrates  with 
special  emphasis  on  the  reduction  of  surface  wave  losses  using  the  finite  difference- time 
domain  technique.  A  microstrip  patch  antenna  is  chosen  as  the  radiating  element  since 
it  is  a  popular  choice  for  a  low  profile  antenna  which  is  used  with  a  substrate  [147]- 
[152],[111],[139]-[144].  In  designing  microstrip  antennas,  a  few  of  the  key  features, 
which  are  commonly  optimized,  are:  minimizing  the  overall  antenna  size,  maximizing 
the  antenna’s  bandwidth  in  frequency,  and  minimizing  the  losses  due  to  surface  waves. 
The  overall  antenna  size  for  operation  at  a  specified  frequency  range  can  typically 
be  reduced  by  using  a  higher  dielectric  constant  substrate.  The  bandwidth  can  be 
increased  by  increasing  the  height  of  the  microstrip  patch.  However,  these  two  changes 
tend  to  increase  the  loss  due  to  surface  waves,  which  is  undesirable. 
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The  presence  of  guided  waves  in  the  substrate  is  undesirable  for  the  following 
reasons.  First  of  all,  power  coupled  into  guided  waves  is  power  which  is  not  radiated 
which  results  in  a  lower  radiation  efficiency  for  the  antenna  [140].  Also  with  a  finite 
size  substrate,  the  interaction  of  the  guided  waves  with  termination  of  the  substrate 
could  cause  unwanted  high  side  lobes  or  back  lobes.  When  phased  arrays  are  mounted 
on  a  common  substrate,  significant  coupling  between  elements  can  occur  which  can 
cause  limitations  in  the  scanning  range  of  the  array  [147], [148]. 

The  uniform  substrate  is  altered  by  removing  a  rectangular  ring  of  the  dielectric 
substrate  which  surrounds  the  radiating  element,  and  is  shown  in  Figure  6.1.  By  re¬ 
moving  rings  of  dielectric  from  the  substrate  which  surround  the  microstrip  patch,  the 
surface  or  guided  waves  are  reflected  at  the  interfaces  and  can  couple  more  efficiently 
into  radiated  modes. 

In  Section  6.2  a  brief  discussion  of  the  finite  difference-time  domain  technique  is 
given.  A  perfectly  matched  layer  absorbing  boundary  condition  is  used  because  of  its 
superior  performance  over  other  absorbing  boundary  conditions  and  is  described  in 
detail  in  Section  6.2.  The  results  are  shown  and  discussed  in  Section  6.3.  In  Section 
6.4  the  usefulness  and  effectiveness  of  the  non-uniform  substrate  in  terms  of  reducing 
the  surface  wave  generation  is  discussed  and  additional  configurations  and  uses  of 
non-uniform  substrates  are  commented  on. 
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Figure  6.1:  Microstrip  patch  antenna  on  a  non-uniform  dielectric  substrate. 
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6.2  Method  of  Solution 

6.2.1  Finite  Difference  —  Time  Domain  Technique 

The  finite  difference-time  domain  (FD-TD)  technique  is  used  to  analyze  this  problem 
[1]-[12].  This  method  was  chosen  because  it  can  be  used  to  analyze  multiple  frequen¬ 
cies  with  a  single  simulation,  and  it  is  a  good  choice  for  low  loss  models  and  rectilinear 
geometries. 

A  rectangular  three-dimensional  grid  was  used  as  the  computational  domain. 
Spatially,  the  computational  domain  is  discretized  such  that  the  largest  grid  dimen¬ 
sion  is  less  than  about  a  twentieth  of  the  shortest  wavelength  of  interest.  The  time 
step  approximately  must  be  less  than  the  smallest  grid  dimension  divided  by  the 
speed  of  light  in  order  to  guarantee  the  stability  of  the  algorithm.  The  difference 
equations  which  are  utilized  to  update  the  electric  and  magnetic  fields  are  obtained 
from  Maxwell’s  equations  in  differential  form  by  approximating  spatial  and  temporal 
derivatives  with  center  differences.  Electric  and  magnetic  fields  are  discretized  on 
Yee’s  lattice  [1].  Initial  conditions,  which  are  imposed  at  the  beginning  of  the  simula¬ 
tion,  are  that  all  electric  and  magnetic  fields  start  at  zero.  At  each  subsequent  time 
step,  fields  are  imposed  at  appropriate  locations  in  order  to  simulate  the  excitation 
source.  The  imposed  electric  field  has  a  time  dependence  of  a  Gaussian  pulse  modu¬ 
lated  by  a  sine  wave,  which  contains  energy  in  a  band  of  frequencies  allowing  multiple 
frequencies  to  be  analyzed  simultaneously.  Boundary  conditions  are  imposed  at  per¬ 
fectly  conducting  surfaces  and  at  the  outer  boundary.  At  the  perfectly  conducting 
surfaces,  the  tangential  electric  field  is  set  to  zero.  Absorbing  boundary  conditions 
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must  be  applied  at  the  outer  boundary  of  the  computational  domain  in  order  to  sim¬ 
ulate  unbounded  space.  Perfectly  matched  layer  absorbing  boundary  conditions  are 
used  at  the  outer  boundary  [78]-[83].  Finally,  by  applying  the  difference  equations 
obtained  from  Maxwell’s  equations  in  addition  to  the  appropriate  initial  and  bound¬ 
ary  conditions,  the  electric  and  magnetic  fields  can  be  calculated  everywhere  within 
the  computational  domain  for  each  time  step. 

The  source  used  with  this  problem  was  a  coaxial  line  probe  exciting  the  mi¬ 
crostrip  patch.  The  treatment  of  the  source  condition  is  discussed  in  Subsection 
4.3.1.  The  equation  governing  the  update  of  the  electric  field  at  the  source  is  given 
in  Equation  4.1.  The  calculation  for  the  radiation  impedance  of  the  antenna  and  the 
current  and  voltage  on  the  coaxial  probe  are  also  discussed  in  Subsection  4.3.1. 

The  calculations  for  total  radiated  and  guided  power  are  done  using  time  domain 
results  from  the  FD-TD  simulation.  First,  the  electric  and  magnetic  fields  are  calcu¬ 
lated  everywhere  in  the  computational  domain  and  then  the  tangential  electric  and 
magnetic  fields  on  a  surface  enclosing  the  antenna  are  recorded.  The  fields  are  then 
Fourier  transformed  to  obtain  frequency  domain  data  and  the  time  average  Poynting 
power,  (5(/)),  is  calculated. 

(s(/))  =  5Re(®(/)  X  iTif)]  (6.1) 

The  total  guided  power  is  then  approximated  by  integrating  the  Poynting  power 
passing  through  the  integration  surface  which  lies  within  the  substrate.  The  inte¬ 
gration  surface  is  placed  far  enough  from  the  source  such  that  most  of  the  power 
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in  the  substrate  is  presumably  guided.  The  total  radiated  power  is  then  calculated 
by  integrating  the  Poynting  power  passing  through  the  integration  surface  above  the 
substrate. 


6.2.2  Perfectly  Matched  Layer  Absorbing  Boundary  Con¬ 
dition 

The  perfectly  matched  layer  (PML)  absorbing  boundary  condition  (ABC)  [78]-[83] 
was  recently  introduced  in  1994  by  Berenger,  and  was  chosen  to  be  used  over  the 
standard  one-way  equation  absorbing  boundary  conditions  (Section  2.11)  [67]-[77] 
because  of  its  superior  performance  in  general  and  especially  with  dielectric  substrates. 

The  basic  premise  of  the  PML  ABC  is  to  implement  a  media  which  is  reflec¬ 
tionless  for  waves  incident  from  a  free  space  or  dielectric  region,  and  which  is  highly 
lossy.  By  providing  sufficient  loss  in  the  PML  layer  the  truncation  of  the  computa¬ 
tional  domain  can  then  be  performed  by  using  perfect  electric  conductor  and  still  not 
introduce  significant  reflections. 

In  order  to  facilitate  the  analysis  of  the  perfectly  matched  layer  absorbing  bound¬ 
ary  condition,  a  two  dimensional  case  will  be  considered  first.  Then  the  two  dimen¬ 
sional  treatment  will  be  extrapolated  to  three  dimensions. 

The  electromagnetic  fields,  in  a  PML  for  a  two  dimensional  case  where  there  is 
no  variation  in  the  2  direction,  are  governed  by  the  following  equations  for  the  TM 
polarization. 
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erCo' 


dR. 


dt 


-j“  (Tj^Ey 


dx 


flrflQ 


dt 

dH,y 

dt 


“I"  ^mxHzx 


(^myHzy 


_d^ 

dx 

dEg; 

dy 


(6.3) 

(6.4) 

(6.5) 


Note  that  there  are  four  field  components  {E^,  Ey,  Hzx  and  Hzy)  instead  of  the 
standard  three.  The  sum  of  Hzx  and  Hzy  can  be  interpreted  as  the  normal  Hz. 
ax  and  ay  can  be  interpreted  as  anisotropic  electric  conductivities,  a^x  and  amy 
have  a  similar  interpretation  as  magnetic  conductivities  with  respect  to  Hzx  and  Hzy\ 
however  they  have  a  less  clear  relationship  with  Hz.  Note  that  when  =  cfmy, 
the  PML  equations  reduce  to  Maxwell’s  equations  and  the  PML  medium  becomes  a 
normal  dielectric  or  magnetic  material. 

The  following  conditions  are  imposed  on  the  properties  of  the  PML  layer  in 
order  to  make  it  reflectionless: 


_  ^mx 

6r6o 


(6.6) 


^y  _  ^my 

CrCo 


(6.7) 


These  conditions  are  essentially  identical  to  those  imposed  for  the  impedance  match¬ 
ing  problem  for  a  wave  at  normal  incidence.  It  will  be  shown  that  with  the  alterations 
to  Maxwell’s  equations  outlined  above  that  reflectionless  transmission  can  be  achieved 
for  all  angles  of  incidence  once  additional  constraints  on  the  PML  parameters  are  im- 
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The  continuity  of  tangential  electric  and  magnetic  fields  at  a  dielectric/PML 
layer  interface  can  be  established  by  integrating  the  sum  of  Equations  6.4  and  6.5  on 
the  surface  Si  shown  in  Figure  6.2.  Applying  Gauss’  theorem,  the  following  equation 
is  obtained: 

J^VxE-dS  =  J^E-dl  =  /Xr2//0  +  amyH.y 

f  dH^y  +  Hzx  (I, 

+  -  (6-8) 

Making  Ay  of  Si  vanishingly  thin,  the  surface  integral  of  Hzy  and  Hz^  vanishes,  and 
it  is  apparent  that 


^dielectric^^  -  El"'^Nx  =  0 


(6.9) 


^dielectric  _ 

Similar  arguments  apply  to  show  that  other  components  of  the  electric  and  magnetic 
fields  which  are  tangential  to  the  interface  are  continuous. 

Since  the  tangential  fields  are  continuous,  phase  matching  must  apply,  and  thus 


^pml 


(6.10) 


the  tangential  wavenumber,  kx^  must  be  continuous 
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Figure  6.2:  Integration  surface  for  interface  between  a  dielectric  medium  and  a  PML 
medium. 
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Ki  =  c  =  Kt  =  0  (6.12) 

for  propagating  plane  wave  incidence,  k'^.  and  A:"  are  the  real  and  imaginary  portions 
of  the  tangential  wave  number  and  the  subscripts  r  and  t  denote  incident,  reflected 
and  transmitted  wavenumbers  respectively. 

In  deriving  a  wave  equation  for  a  PML  medium,  the  time-harmonic  equations 
in  a  PML  layer  will  be  used  and  a  time  dependence  is  assumed. 


Ex{-iu;€reo  -1-  Gy) 

=  ikyi^Hzx  “1"  Hzy) 

(6.13) 

Ey{-iujereo  +  Gx) 

=  ikx{Hzx  “1”  Hzy) 

(6.14) 

Hzxi,  “b  ^mx) 

=  ikxi^Hyx  "f"  Hyz) 

(6.15) 

Hzyi-iuHrlJ-O  +  (^my) 

”  Hxz) 

(6.16) 

By  simple  manipulation  of  the  above  equations  a  wave  equation  in  a  PML  medium 
with  Cr,  /jLr,  cTj,  ffy,  (Traxi  <^my  cun  be  derived. 

^2'  (  -|-  Gy  \  ^  ^2  (  -|“  Gx  \ 

y  “b  Gfnx  j  ^  \  ~b  ^my  j 

-t-(-ia;€reo  +  Gy){-i(jjereQ  -t-  ctj)  =  0  (6-17) 


k 


2 

X 


1  -  Gyl{iijjereQ) 
1  ^mx! 


+  kl 


1  -  (Jx! 

1  (^my/ii^l^rPoi) 


(6.18) 


-t3[l  -  c’»/(»w£,eo)][l  -  <Ti/(ia>ereo)l  =  0 
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k 


y 


Cry/(2CJ6r6o)]  [1 


p  \l-(Ty/{iujereo) 
[1-  ax/iiivereo) 


1  ^my  /  l^rl^o) 

1  ~  ^mx/ 


(6.19) 


where  kr  =  uj^/eiT^JI^. 

In  calculating  the  reflection  coefficient  for  a  plane  wave  incident  from  a  dielectric 
region  on  a  PML  layer,  a  plane  wave  with  Hmc  =  is  used.  Using 

Equations  6.11-6.14,  the  equations  in  the  dielectric  region  characterized  by  Cri  and 


-h  H,y)  =  +  Re-^'‘^yy)e^'^'^ 


(6.20) 


E,= 


ih. 


ly 


-lue^eri 


Hoie^'^'yy  - 


(6.21) 


Ey  =  Hoie^'^'yy  + 

-ZU/eoCrl 


(6.22) 


where  R  is  reffection  coefficient  between  the  dielectric  and  PML  regions.  The  equa¬ 
tions  in  the  PML  region  characterized  by  Crz,  A‘r2)  (’’d  ctj,,  amx  and  cr^j,  are 


-F  H,y)  =  THoe^^^^^e^'^^yy  (6.23) 

- J'H^^iki.x^iktyy 

-iuje^er^  +  Oy 


(6.24) 
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Ey  — 


—iki. 


-iu>eoer2  +  o'; 


-ri7oe*'=  tx 


(6.25) 


where  T  is  the  transmission  coefficient  between  the  dielectric  and  PML  regions. 

Setting  the  tangential  electric  and  magnetic  fields  continuous  at  the  boundary 
results  in 


l  +  i2  =  T 


(6.26) 


-R)  =  - T 

—lUJ€QCrl  —itUCQCr^  +  Oy 

Combining  Equations  6.26  and  6.27,  the  reflection  coefficient  is  found  to  be 


(6.27) 


where 


R^\-^ 

1  +p 


(6.28) 


kfy  1 

kiy  ^y! (^^^O^rl) 


(6.29) 


It  is  apparent  that  the  interface  will  be  reflectionless  when  p  =  1. 

Applying  the  conditions  for  a  matched  medium,  Equations  6.6  and  6.7,  the 
normal  transmitted  wave  number  is 


kty  — 


i 


^r2[l  - 


1  -  Oy/ (ia;€oer2) 

1  -  (r^/{iu}eQer2)\ 


(6.30) 
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where  kr2  —  a;^er2€oMr2A^o-  After  some  consideration,  it  can  be  observed  that  reflec¬ 
tionless  transmission  for  all  angles  will  occur  only  under  the  following  conditions: 

Cri  =  er2  (6.31) 

Hri  =  Hr2  (6.32) 

Ox  =  Omx  =  0  (6.33) 


Under  these  conditions, 


kty  =  [1  -  Oyl{iujeQer2)\\lk^i  -  k^x  (6-34) 

and  it  is  observed  that  p  =  1  in  this  case.  To  summarize  the  conditions  for  reflection¬ 
less  transmission,  all  the  parameters  in  the  PML  layer  must  be  identical  to  those  in 
the  dielectric  layer  except  the  loss  factors,  a  and  am  in  the  direction  normal  to  the 
interface.  The  same  conditions  hold  for  interfaces  between  two  PML  layers. 

As  mentioned  earlier,  the  operating  principle  of  the  PML  ABC  is  that  the  PML 
layer  is  reflectionless  and  provides  significant  loss  as  the  wave  propagates  through 
it.  Since  the  PML  layer  is  typically  backed  by  a  perfect  conductor,  the  total  loss 
as  a  wave  propagates  through  the  PML  layer  will  be  determined  by  its  thickness 
and  its  conductivities.  In  order  to  achieve  low  reflections  from  the  conductor  backed 
PML  layer,  large  values  of  conductivity  are  used.  Applying  center  differences  to  the 
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PML  equations,  Equations  6. 2-6. 5,  is  impractical  since  the  grid  cell  would  need  to  be 
greatly  reduced  in  order  to  accommodate  the  desired  high  conductivities.  Instead, 
explicit  exponentially  differenced  equations  are  used  in  order  to  maintain  accuracy 
and  stability. 

The  explicit  exponentially  differenced  equations  are  obtained  from  solving  the 
quasi-static  charge  relaxation  problem.  The  following  equation  for  Ey  is  solved  using 
quasi-static  approximations. 


dE, 


dEy  d{H,,  -h  H,y) 

-I-  — - - - 


(6.35) 


Assuming  that  the  time  constant  associated  with  the  conductivity  is  small  compared 
to  the  period  of  the  frequency  of  interest,  the  magnetic  field  terms  can  be  considered 
slowly  varying  in  time  and  is  approximated  as  constant.  The  homogeneous  solution 
to  Equation  6.35  is 


Eyh{t)  =  Ey{t  = 


(6.36) 


The  particular  solution,  approximating  the  magnetic  field  as  constant  in  time,  is 


p  (f\  _  1  9{Hzx  +  Hzy) 

(yOC 


(6.37) 


Combining  the  homogeneous  and  particular  solutions  results  in 


Ey{t)  =  Ey{t  = 

0  X  dx 


(6.38) 
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By  taking  the  solution  at  time  step  n  + 1  and  subtracting  the  solution  at  time  step  n 
times  the  following  difference  equation  is  obtained. 

+  1/2)  =  +  +  (6.39) 

X  +  1/2,  j  +  1/2)  -  -  1/2,  j  +  1/2)) 

In  theory  the  loss  factor  can  be  chosen  to  be  arbitrarily  large  and  an  arbitrarily 
thin  PML  layer  could  be  used.  However,  in  implementing  the  equations  in  the  FD- 
TD  method,  the  discretized  nature  of  space  must  be  considered  since  the  electric  and 
magnetic  fields  are  separated  by  half  grid  cells.  Hence,  it  is  desirable  to  chose  PML 
layers  which  are  4  to  16  cells  thick,  and  to  grade  the  conductivity  smoothly  from  zero 
to  some  maximum  such  that  discretization  errors  will  be  minimized.  One  manner  of 
grading  the  conductivity  is  by  using  a  polynomial  curve. 

=  (^max[p/(Ap)]”  (6.40) 

where  Ap  is  the  PML  thickness,  p  is  the  position  in  the  PML  layer,  and  n  is  a 
smoothness  factor.  Quadratic  grading  of  the  conductivity  has  shown  itself  to  work 
well.  Note  that  since  the  electric  and  magnetic  fields  are  separated  by  a  half  cell, 
this  half  cell  should  be  accounted  for  in  computing  the  desired  conductivities  for  the 
affected  electric  and  magnetic  fields. 

The  effective  reflection  coefficient  for  a  PML  layer  backed  by  a  perfect  conductor 
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Ox  ^0 


Figure  6.3:  Selection  of  conductivities  for  a  PML  ABC  in  two-dimensions. 


will  be 


where  6  is  the  angle  of  incidence,  9  =  0  degrees  for  normal  incidence,  and  Cr  = 
l/V^rfo/^rMo-  The  reflection  coefficient  is  typically  chosen  to  be  around  i?(0)  =  10“® 
as  a  compromise  between  minimizing  reflections  and  avoiding  numerical  noise.  Once 
the  cell  size,  PML  layer  thickness,  permittivity,  permeability  and  desired  reflection 
coefficient  are  determined,  then  the  conductivities  in  the  PML  layer  can  be  determined 
using  Equations  6.6-6. 7,  and  6.40-6.41.  Figure  6.3  shows  the  selection  of  the  various 
conductivities  needed  to  implement  a  PML  ABC. 

In  extending  the  use  of  the  PML  ABC  to  three-dimensions,  only  a  few  modifi¬ 
cations  are  needed.  All  six  field  components  in  rectangular  coordinates  are  split  in 
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two,  and  a  separate  PML  conductivity  is  used  in  each  direction.  The  resulting  set  of 
twelve  equations  are 


^r^o-^  +  ayE^y 

dExz  „ 


^■^yx  .  jT, 

oj  '  ^x^yx 
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I  (Jz^yz 


\yjL^ZX  1— I 
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^reo-^  +  (TyEzy 
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o.  +  ^my-tixy 


^  2/2;  .  y  Y 

/^r/^0  'q^  t  O'mz-^x 
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l^rt^O  “h  (^mxNy; 


A*rA*0  d"  (^mzHyz 
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d{Hy,  +  Hyz) 
dz 

d{Hz,+Hzy) 

dx 

d{H,y+H,,) 

dz 

d{Hy,  +  Hy,) 
dx 

d{Hz,y-\-H,,) 

dy 

9{Ezx  +  Ezy) 
dy 

9{Dyx  d*  Dyz) 
dz 

d(Ezx  +  Ezy) 
dx 

d{Exy  d~  Exz') 

dz 

d{Eyx  d-  Eyz) 

dx 

d{Exy  +  Exz) 

dy 


(6.42) 


(6.43) 


(6.44) 


(6.45) 


(6.46) 


(6.47) 


(6.48) 


(6.49) 


(6.50) 


(6.51) 


(6.52) 


(6.53) 


These  twelve  equations  for  a  PML  medium  reduce  to  Maxwell’s  equations  when  ax  = 


Gy  a^  and  a^nx  —  —  ^mz- 
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In  selecting  the  PML  layer  parameters  the  same  rules  are  used  for  three-dimen¬ 
sions  as  for  two-dimensions.  The  summary  of  the  rule  is  that:  all  parameters  in  the 
PML  layer  must  be  identical  to  the  dielectric  layer  except  for  the  conductivities  in  the 
direction  normal  to  the  interface.  This  same  rule  still  applies  to  interfaces  between 
two  PML  layers. 

In  setting  up  the  PML  ABC  surrounding  a  free  space  region  in  three-dimensions, 
the  selection  of  conductivities  will  be  as  follows:  (1)  the  faces  of  the  computational 
domain  will  have  loss  only  in  the  direction  normal  to  them,  (2)  the  edges  will  have 
loss  in  the  two  directions  normal  to  them,  and  (3)  the  corners  will  have  loss  in  all 
three  directions. 

The  parameters  for  a  PML  layer  which  is  adjacent  at  a  constant  x  plane  to 
a  dielectric  medium  with  parameters  Cr,  fir,  cr  =  0,  and  =  0,  are  Cr,  fir,  Oy  = 
O,  =  Cmy  =  (Jmz  =  0,  and  ax/(ereo)  =  (TmxlifJ^rPo)  =  cTo/(erCo)  #  0.  The  difference 
equations  in  the  PML  layer  at  the  -\-x  face  of  the  computational  domain  are 


E^y\i  +  hh  =  ^xy{i  +  h  h  k)  +  m 


At 


erAy 


Hzt^i  +  lj  +  lk) 


(i  +  i  i  it)  -  -ff- ‘  (i  +  -  5.  +  5-  ^  i 


(6.54) 


(6.55) 
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'Ho 


y,rAz 


(6.61) 


X  +  1)  +  £;.(«, i  +  1. 1  +  1)  -  j  +  1,  *)  -  £J.(i,  J  +  1.  A:)] 


+  \,j,  A  +  1)  =  +  1,1,  fc  +  l)e-'’«'i‘/"‘"  - 

X  + 1,1,  +  5)  +  E^y{i  +  1,1,  A  +  1) 


€oer(e-‘"o^‘/‘-^o  -  1) 


GaHrUdAx 


-E;,(i,j,k  +  k)-E:^(ij,k  +  k)] 


(6.62) 
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Hyz  *  (z  +  j,  A:  +  i)  =  Hyz  2 +  i  /j  +  i)  _  J_ 


% 


Ar 


firAkZ 


(6.63) 


^  [^xyi^  +  2>-^>  ^  +  1)  +  E^zi'^  +  +  I)  —  E^y{i  +  k)  —  E^^(i  +  I,j,  A;)j 

+  i,  J  +  i,  k)  =  Hzxhi  +  hj  +  h  + 

X  [Eyxi'i'  +  l,i  +  k)  +  Ey^{i  +  1,  j  k) 


€oer(e-‘^»^*/^’-*<>  -  1) 
aofXrHo^X 


(6.64) 


-E;xihj  +  hk)-E^,{i,j  +  lk)] 


Hzy^  (^’  +  |,  i  +  5,  k)  =  Hzy  +  2  >  •?  +  2>  ^)  +  ~ 


Vo 


At 


2fj,rAy 


(6.65) 


X  K(>  +  |,  j  +  1,  *:)  +  +  |,  J  +  l,k)  -  +  1,;,  i)  -  +  1.  j,  *:)] 


Note  that  these  equations  for  the  split  field  quantities  can  be  used  for  non-PML 
media,  and  the  standard  field  quantities  can  be  obtained  simply  by  adding  the  split 
quantities. 


6.2.3  Radiation  Pattern  over  a  Multi-Layer  Media 

In  order  to  calculate  the  far-field  radiation  pattern,  it  is  necessary  to  use  a  multi-layer 
Green’s  function  [106], [153].  For  this  configuration,  there  are  three  layers  with  the 
top  layer  being  free  space,  the  middle  layer  being  a  dielectric  with  €r,  and  the  lower 
layer  being  a  perfect  conductor  as  shown  in  Figure  6.4.  The  source  is  located  in  the 
upper  two  layers  and  the  far-field  electric  fields  are  only  calculated  in  the  free  space 
region.  The  mathematical  formulation  of  Huygens’  principle  for  this  configuration. 
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Figure  6.4:  Three  layer  configuration.  A  dielectric  substrate  on  a  perfect  conductor. 


assuming  a  e  time  dependence,  for  the  electric  field,  is 

•  h  X  H(f')  +  V  x  ^oo^(f,r^)  •  h  x 
+  J Ig,  •  n  X  H{f')  +  V  X  •  n  x  E(r')| 

(6.66) 

where  S[  +  5^  form  a  closed  surface  with  the  ground  plane,  S'l  is  the  portion  in  the 
free  space  region,  and  S2  is  the  portion  in  the  dielectric  layer.  The  subscripts  to 
the  Green’s  function.  Gab,  are  a  for  the  observation  layer,  and  b  for  the  source  layer. 
The  superscripts,  ef  f  and  mff,  to  the  Green’s  function  indicate  that  the  Green’s 
function  is  only  valid  in  the  far-field  and  that  it  is  applied  to  the  electric  field  and  the 
magnetic  field  respectively,  f  is  the  far-field  observation  location  and  f'  is  the  source 
location.  Note  that  the  observation  location  is  always  higher  in  z  than  the  source 


location. 
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The  Green’s  function,  ,  in  the  spectral  domain,  is 


Goo  {r>,f'<)  =  ^  r  dK  r 

OTT  7—00  7—00  Kqz 

{ +  O(^0,)A(-4z)e“”''«"'']  (6.67) 


+ 


'v{koz)v{koz)e^'^°‘^’‘-^''>  +  Rl^v{koz)v{-ko,)e^'^°^^^+^'^] } 


where  the  wavenumbers  are  defined  as 


kj  -  (jJy/CjHj 


(6.68) 


and  j  is  the  number  of  the  layer,  where  0  is  the  free  space  region,  1  is  the  dielectric 
region,  and  2  is  the  perfect  conductor. 


kj  =  xkx  +  yky  +  zkjz  (6.69) 

=  \ll=]  -kl-kl  (6.70) 

p  =  xx  +  yy  (6-71) 

kp  =  xkx  +  yky  (6.72) 
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The  horizontal  and  vertical  unit  vectors  are  defined  as 


kj  X  z 
kj  X  z 


=  -  y^x) 

Kp 


"t"  y^y)  I  ^ 


kjkp 


h{  kjz)  —  h(^kjz) 


f./  \  _  ^jzi&kx  +  yky)  ^kp 


The  reflection  coeflScient  looking  down  from  layer  i  is 


The  reflection  coefficient  looking  up  from  layer  (.  is 


where  s  represents  either  the  TE  or  TM  polarization.  The  reflection 
looking  from  layer  i  to  layer  j  are 


^TE  _  l-^j^iz  y’jkjz 

I^jkiz  ”1“  y'ikjz 


(6.73) 

(6.74) 

(6.75) 

(6.76) 

(6.77) 

(6.78) 

coefficients 


(6.79) 


• 
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pTM  ^j^iz  ^ikjz 

(6.80) 

• 

z  =  r  cos  $ 

(6.81) 

• 

p  —  rsmO 

(6.82) 

• 

X  =  pCOS(f) 

(6.83) 

• 

y  =  psin0 

(6.84) 

• 

hz  =  y/ko  -  =  ko  COS7 

(6.85) 

• 

kp  =  y/k^  -  k^^  =  kosiwy 

(6.86) 

kx  =  kp  cos  5 

(6.87) 

• 

ky  =  kp  sin  J 

(6.88) 

1 

• 

• 
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Changing  the  integrals  of  Equation  6.67  into  cylindrical  coordinates  gives 


r  ^Se^k,p[cosSco.4>+B\nSsin<l>]^iy^ 

OTT  «/  0  »/  0 


j_ 

koz 


{ +  RjSHh.)H-h,)e"‘'“‘’]  (6.89) 


Approximating  the  above  dS  integral  using  the  stationary  phase  approximation  at 
S  =  4>  results  in 


(6.90) 

87r2  Jo  ^  V  Kp  ^oz 

{[h(/i:o2)h(/i:o2)e“**'‘>"^'  +  R^oh{ko^)h{-koz)e'''°^'''] 

+  [v{koz)v{koz)e~"''°‘’''  +  i2™0(A;o2)?)(-A:oz)e**'“"*'] 


Approximating  the  above  dkp  integral  using  the  stationary  phase  approximation  at 
kp  =  ko  sin  9  results  in 

3:"(.>y<)  =  (e.91) 

{ [h{ko.)h{h.)e-^-’‘  +  It^hih,)h(-k,,)e'«”'] 

L  ^  d— “(pj/Cp  — ruQ  Sin,  0 
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And  simplifying,  results  in 


^// 


Goo 


f>ikor  _ 

_ _ Q-ikp-p' 

Airr 


(6.92) 


{ [h{koMko,)e-^’‘<^^^'  +  R'^fh{koM-ko,)e^’^^^^'] 

'■  ^  J  d=<l>,kp=koSin 


With  the  far-field  approximations,  the  solution  is  only  valid  for  observations  points 
in  the  far-field  in  layer  0  and  not  in  the  near-field  of  layer  1.  The  associated  Green’s 
function  for  a  source  in  layer  1  and  a  far-field  observation  point  in  layer  0  is 


G..  <'■>.’■<)  = 


(6.93) 


VTE 

^Ul^O 


+ 


1  - 

ytm 

-Aj  n  , 


h{koz)h{koz)e-^’^°^^'  +  Rj,fe'^^^^‘^^h{koz)hi-koz)e^^°^^'] 


where  the  upward  transmission  coefficients  from  layer  m  to  layer  £  with  m  >  £  are 
defined  as 


yTE 

^Um— 


YTE 


(l  +  iZJc^+i)) 

(1  +  i2TEe2ifc/zA«) 


(6.94) 


y-TM  _ 


ke  k(^e-\-i)z 
kiz  ^(£+1) 


Y™ 


(1  -  e2**=^^'u) 


(6.95) 
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The  upward  transmission  coefficients  for  adjacent  layers  are  defined  as 


X, 


TE 

Um— >m— 1 


( 


1  +  ^S)) 


(6.96) 


km-1  kmz _ (l  -^Um) _ 


(6.97) 


The  Green’s  function  applied  to  the  electric  field  is  identical  to  the  one  for  the 
magnetic  field,  except  that  the  horizontal  direction  is  now  associated  with  the  TM 
polarization  and  the  vertical  direction  with  the  TE  polarization. 

=eff 

Goo{f>,f<)  =  (6-98) 

{[h{koz)h{koz)e~'''°^'‘'  +  I^Qh{koz)h{-ko:z)e’‘^°"''’] 

+  [v{koz)v{kQz)e~^''’'"^’  +  R™v{koz)v{-koz)e'^°""']} 


The  associated  Green’s  function  for  a  source  in  layer  1  and  a  far-field  observation 
point  in  layer  0  is 


^ffr-  _,v  koz 


^or(r>y<)  =  ^e-V^ 


(6.99) 


{yTM 

v-TE 
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6.3  Results  and  Discussion 

In  performing  the  FD-TD  simulations  on  the  microstrip  antenna  and  the  vertical 
and  horizontal  dipole  configurations,  the  grid  cell  which  was  used  had  dimensions  of 
Arc  =  Ay  =  1.25  cm,  and  Az  =  1.2  cm.  The  PML  ABC  boundary  condition  was 
chosen  to  be  6  cells  thick  in  all  directions,  a  quadratic  grading  of  the  conductivity 
was  used,  and  a  reflection  coefficient  of  10”^  was  chosen  for  the  PML  ABC  at  normal 
incidence. 

Before  examining  a  microstrip  patch  antenna  configuration,  two  simple  dipole 
cases  are  examined.  The  first  is  a  vertical  dipole  configuration  in  a  non-uniform 
substrate  as  shown  in  Figure  6.5  which  is  compared  with  an  identical  vertical  dipole 
in  a  uniform  substrate. 

Three  vertical  dipole  configurations  are  examined,  and  all  three  have  the  same 
dipole  height  and  substrate  characteristics.  One  configuration  is  for  a  uniform  sub¬ 
strate  and  the  other  two  are  for  non-uniform  substrates,  h  is  the  thickness  of  the 
substrate  and  the  height  of  the  vertical  dipole. 


h  =  4.8  cm 


(6.100) 


The  dielectric  constant  of  the  substrate  is 


e,.  =  12.8 


(6.101) 


Both  non-uniform  substrate  configurations  have  the  vertical  dipole  located  at 
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Figure  6.5:  Vertical  dipole  in  a  non-uniform  dielectric  substrate. 
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the  center  of  the  center  block  of  dielectric  substrate.  Wcxi)  Wcyi,  Wcx2  '^cy2  ^re 
the  dimensions  of  the  center  block  of  dielectric  substrate  for  the  first  and  second 
non-uniform  substrate  configurations. 

Wcxi  =  t^cyi  =  10  cm  (6.102) 

Wcx2  -  Wcy2  =  7-5  cm  (6.103) 

Wrxu  Wryi,  Wrx2  and  Wry2  are  the  widths  of  the  free  space  ring  surrounding  the  center 
block  of  dielectric  for  both  configurations. 

Wrxi  =  =  17.5  cm  (6.104) 


Wri2  =  Wry2  =  11-25  Cm 


(6.105) 


These  widths  were  chosen  to  be  between  a  quarter  and  half  a  wavelength  in  size  in 
order  to  allow  the  power  to  radiate  into  the  free  space  region. 

The  surfaces  waves  are  guided  modes  and  have  cut-off  frequencies  of 


(6.106) 


where  n  =  0, 1, 2, 3, ...,  with  n  even  for  TM  modes  and  n  odd  for  TE  modes  [106], [140]. 
Co  is  the  speed  of  light,  h  is  the  substrate  height,  the  relative  permittivity  and  Hr 
the  relative  permeability  of  the  substrate.  So  the  TMq  mode  will  always  be  present 
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since  it  is  never  cut-off,  and  for  this  case,  the  cut-off  frequency  for  the  TEi  mode  is 
455  MHz,  and  for  the  TM2  mode  is  910  MHz. 

Figure  6.6  shows  the  radiation  impedance  of  all  three  of  the  vertical  dipole  con¬ 
figurations.  The  peaks  in  the  radiation  impedance  occur  at  940  MHz  for  the  uniform 
substrate  case,  at  835  MHz  for  the  first  non-uniform  substrate  case,  and  at  935  MHz 
for  the  second  non-uniform  substrate  case.  It  is  noted  that  2  guided  wave  modes  can 
propagate  in  the  first  non-uniform  substrate  case  and  3  guided  modes  can  propagate 
in  the  other  two  cases.  It  is  observed  from  this  that  the  resonant  frequency  can  be 
shifted  by  using  different  size  free  space  rings  surrounding  the  antenna.  These  reso¬ 
nances  correspond  to  the  full  wavelength  resonance  of  a  dipole,  or  a  half-wavelength 
resonance  of  a  monopole  over  a  ground  plane.  The  corresponding  bandwidths  are  260 
MHz,  90  MHz,  and  95  MHz  for  the  uniform,  first  and  second  non-uniform  substrate 
cases.  The  addition  of  the  non-uniformity  has  significantly  increased  the  configura¬ 
tion’s  Q  which  results  in  much  lower  bandwidths. 

The  approximate  percentage  of  power  guided  by  these  configurations  are  shown 
in  Figure  6.7.  The  surface  on  which  the  Poynting  power  is  integrated,  is  always  chosen 
to  enclose  the  antenna  and  the  free  space  ring.  The  surface  is  approximately  2-4  cells 
larger  in  the  x  and  y  directions  than  the  free  space  ring  for  all  the  cases  considered. 
The  approximate  percentages  of  power  guided  by  these  three  configurations  over 
the  frequency  range  surrounding  their  resonant  frequencies  are  98%  for  the  uniform 
substrate  case,  22%  for  the  first  non-uniform  substrate  case,  and  30%  for  the  second 
non-uniform  substrate  case.  The  reduction  with  both  non-uniform  cases  is  significant 
compared  to  the  uniform  case  over  the  frequency  range  around  the  resonant  frequency. 
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Figure  6.6:  Radiation  impedance  of  a  vertical  dipole  in  a  uniform  and  two  non-uniform 
substrate  configurations. 
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Figure  6.7:  Guided  power  of  a  vertical  dipole  in  a  uniform  and  two  non-uniform 
substrate  configurations. 


The  gain  patterns  for  the  three  configurations  at  their  respective  resonant  fre¬ 
quencies  in  the  x-z  plane  are  shown  in  Figure  6.8.  Only  the  vertical  polarization  is 
plotted  since  the  gain  for  the  horizontal  polarization  is  negligible.  There  is  a  difference 
in  the  overall  levels  of  the  three  gain  patterns  since  significantly  different  amounts  of 
power  are  trapped  in  the  guided  wave  modes  for  each  case.  The  highest  gain  pattern 
is  for  the  first  non-uniform  case  which  has  the  lowest  amount  of  guided  waves  and  the 
lowest  gain  pattern  is  for  the  uniform  case  which  has  almost  all  of  its  power  trapped 
in  the  substrate. 


GAIN 
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Figure  6.8:  Gain  pattern  for  a  vertical  dipole  in  a  uniform  and  two  non-uniform 
substrate  configurations  in  the  x-z  plane  for  the  vertical  polarization. 
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The  second  dipole  case  examined  is  a  horizontal  dipole  on  a  non-uniform  sub¬ 
strate  and  is  shown  in  Figure  6.9.  The  dimensions  of  this  configuration  are  as  follows: 
h  is  the  thickness  of  the  substrate. 


h  =  4.8  cm 

The  dielectric  constant  of  the  substrate  is 


(6.107) 


Cr  =  12.8 


(6.108) 


The  horizontal  dipole  lies  on  top  of  the  substrate  and  is  centered  in  the  dielectric 
box.  Wfid  is  the  width  of  the  dipole. 


Whd  —  15  cm  (6.109) 

Wcx  and  Wcy  are  the  dimensions  of  the  center  block  of  dielectric  substrate. 


Wcx  =  25  cm 


(6.110) 


Wcy  =  10  cm 


(6.111) 


Wrx  and  Wry  are  the  widths  of  the  free  space  ring  surrounding  the  center  block  of 
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Figure  6.9:  Horizontal  dipole  on  a  non-uniform  dielectric  substrate. 
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dielectric. 


Wrx  =  Wry  =  15  ClTl  (6.112) 

Figure  6.10  shows  the  radiation  impedance  of  the  horizontal  dipole  over  uniform 
and  non-uniform  substrates.  For  the  uniform  substrate  case,  the  first  peak  in  radiation 
impedance  occurs  at  437  MHz,  while  for  the  non-uniform  substrate  the  first  peak  is 
at  447  MHz.  These  resonances  correspond  to  a  full  wavelength  resonance.  The 
bandwidth  of  the  uniform  case  is  198  MHz  and  the  bandwidth  for  the  non-uniform 
case  is  25  MHz.  Again  the  addition  of  the  non-uniformity  is  observed  to  significantly 
reduce  the  antenna’s  bandwidth. 

The  approximate  percentages  of  power  guided  by  these  configurations  are  shown 
in  Figure  6.11.  The  approximate  percentages  of  power  guided  by  these  configurations 
at  their  resonant  frequencies  are  41%  for  the  uniform  substrate  case,  and  15%  for 
the  non-uniform  substrate  case.  The  reduction  with  the  non-uniform  case  is  signif¬ 
icant  compared  to  the  uniform  case  over  the  frequency  range  around  their  resonant 
frequencies. 

The  gain  patterns  for  the  horizontal  dipole  configurations  at  their  respective 
resonant  frequencies  in  the  x-z  plane  are  shown  in  Figure  6.12  and  the  gain  patterns 
in  the  y-z  plane  are  shown  in  Figure  6.13.  Only  the  vertical  polarization  is  plotted  in 
the  x-z  plane  since  the  horizontal  polarization  is  negligible,  and  only  the  horizontal 
polarization  is  plotted  in  the  y-z  plane  since  the  vertical  polarization  is  negligible. 
The  gain  for  the  non-uniform  substrate  case  is  about  twice  the  uniform  substrate  case 
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Figure  6.10:  Radiation  impedance  of  a  horizontal  dipole  on  a  uniform  and  on  a  non- 
uniform  substrate. 
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Figure  6.11:  Guided  power  of  a  horizontal  dipole  on  a  uniform  and  on  a  non-uniform 
substrate. 
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Figure  6.12:  Gain  pattern  for  a  horizontal  dipole  on  a  uniform  and  on  a  non-uniform 
substrate  in  the  x-z  plane  for  the  vertical  polarization. 
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Figure  6.13:  Gain  pattern  for  a  horizontal  dipole  on  a  uniform  and  on  a  non-uniform 
substrate  in  the  y-z  plane  for  the  horizontal  polarization. 
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which  is  consistent  with  the  guided  power  associated  with  each  case. 

Having  shown  that  the  non-uniform  substrate  can  reduce  losses  due  to  surface 
waves  with  simple  dipoles,  a  microstrip  patch  antenna  with  an  identical  substrate  is 
analyzed.  The  investigated  microstrip  patch  antenna  over  a  non-uniform  substrate 
is  shown  in  Figure  6.1  and  has  the  following  dimensions,  h  is  the  thickness  of  the 
substrate  and  the  height  of  the  microstrip  patch. 

h  =  4.8  cm  (6.113) 

The  dielectric  constant  of  the  substrate  is 

Cr  =  12.8  (6.114) 

Wax  and  Way  are  the  dimensions  of  the  microstrip  patch. 

Wax  =  10.0  cm  (6.115) 

Way  =  12.5  cm  (6.116) 

Wpx  and  Wpy  give  the  location  of  the  probe  excitation  from  the  center  of  the  microstrip 
patch. 
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Wpy  =  0  cm  (6.118) 

Wcx  and  Wcy  are  the  dimensions  of  the  center  block  of  dielectric  substrate. 

Wcx  =  20  cm  (6.119) 

Wcy  =  22.5  cm  (6.120) 

Wrx  and  Wry  are  the  widths  of  the  free  space  ring  surrounding  the  microstrip  patch 
antenna. 


Wrx  =  20  cm  (6.121) 

Wry  =  20  cm  (6.122) 

The  radiation  characteristics  of  this  microstrip  antenna  on  a  non-uniform  sub¬ 
strate  are  compared  with  an  identical  microstrip  antenna  on  a  uniform  substrate. 
Figure  6.14  shows  the  radiation  impedance  of  the  two  configurations.  The  resonances 
for  the  antenna  on  the  non-uniform  substrate  are  slightly  higher  since  there  is  more 
free  space  in  close  proximity  to  the  antenna.  The  resonant  frequencies  are  determined 
by  the  dimensions  of  the  antenna  and  the  effective  permittivity  surrounding  it,  which 
in  turn,  is  some  average  of  the  dielectrics  surrounding  the  antenna.  The  higher  the 
permittivity,  the  denser  the  medium  is  electrically,  which  results  in  a  lower  resonant 
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frequency  for  this  antenna  configuration.  Since  the  antenna  over  the  uniform  sub¬ 
strate  will  have  a  higher  effective  permittivity  than  the  antenna  over  the  non-uniform 
substrate,  it  has  slightly  lower  resonant  frequencies. 

From  Figure  6.14,  it  is  observed  that  the  first  two  resonances  of  the  antenna  over 
the  uniform  substrate  occur  at  295  MHz  and  at  409  MHz,  and  that  the  corresponding 
resonances  for  the  antenna  over  the  non-uniform  substrate  occur  at  300  MHz  and  421 
MHz.  At  these  frequencies  only  the  lowest  order  surface  mode,  the  TMq  mode,  will 
propagate.  While  the  resonant  frequencies  were  only  slightly  altered,  the  bandwidth 
was  significantly  changed.  The  bandwidth  for  the  first  resonance  changed  from  15 
MHz  to  9  MHz  with  the  addition  of  the  free  space  ring,  and  the  second  resonance 
changed  from  35  MHz  to  11  MHz.  This  reduction  in  bandwidth  is  consistent  with 
the  fact  that  the  free  space  ring  forms  an  additional  dielectric  box  resonator  which 
tends  to  increase  the  structure’s  Q  which  is  inversely  proportional  to  an  antenna’s 
bandwidth. 

The  first  resonance  corresponds  to  the  microstrip  patch  resonating  in  the  x 
direction.  The  first  patch  resonance  occurs  when  the  patch  is  approximately  half  a 
wavelength  long  in  the  substrate  in  the  x  direction.  This  resonance  in  the  x  direction 
occurs  at  a  lower  frequency  than  the  first  one  in  the  y  direction  because  the  patch 
is  excited  by  a  probe  which  is  only  off-center  in  the  x  direction.  Hence,  the  lowest 
resonance  in  the  x  direction  occurs  at  half  a  wavelength  (anti-symmetric  mode),  while 
the  lowest  resonance  in  the  y  direction  occurs  at  a  full  wavelength  (symmetric  mode). 
The  magnetic-wall  approximation  without  considering  fringing  fields  (see  Section  3.4, 
Equation  3.3  [116], [111])  places  the  first  resonance  at  419  MHz.  However,  the  effect 
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Figure  6.14:  Radiation  impedance  of  a  microstrip  patch  antenna  over  a  uniform  and 
over  a  non-uniform  substrate. 
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of  the  fringing  fields,  which  in  this  case  is  quite  significant  due  to  the  large  height  of 
the  patch,  will  reduce  the  first  resonance  to  the  observed  295  MHz  for  the  uniform 
substrate  case  and  300  MHz  for  the  non-uniform  substrate  case. 

The  second  resonance  corresponds  to  a  vertical  resonance,  where  the  patch  is 
approximately  a  half-wavelength  in  height.  The  resonance  is  lower  than  the  resonance 
associated  with  a  wire  without  a  patch  since  the  capacitive  loading  of  the  patch 
reduces  the  resonant  frequency  in  this  case. 

The  amounts  of  power  coupled  into  the  substrate  or  surface  waves  for  the  uni¬ 
form  and  non-uniform  substrate  cases  are  shown  in  Figure  6.15.  The  percentage  of 
guided  power  for  the  uniform  substrate  case  at  its  resonant  frequency  of  295  MHz 
is  3.9%  and  for  the  non-uniform  substrate  at  its  resonant  frequency  of  300  MHz 
is  2.7%,  which  is  approximately  a  30%  reduction  in  the  amount  of  power  coupled 
into  guided  modes.  At  slightly  higher  frequencies,  the  percentage  of  guided  power 
greatly  increases  for  both  configurations,  however  the  non-uniform  substrate  case  has 
markedly  less  guided  power  than  the  uniform  substrate  case.  For  example,  the  second 
resonance  of  the  uniform  substrate  case  of  409  MHZ  has  39%  of  the  power  guided 
while  the  second  resonance  of  the  non-uniform  substrate  case  at  421  MHZ  has  only 
,  14%  of  the  power  guided.  Upon  examining  Figure  6.15,  it  is  observed  that  inserting 
the  discontinuity  into  the  substrate  can  lower  the  loss  due  to  surface  waves  over  a 
wide  range  of  frequencies,  at  least  from  275  MHz  to  500  MHz  for  this  case. 

The  gain  for  the  vertical  polarization  in  the  x-z  plane  is  plotted  in  Figure  6.16 
for  both  configurations  for  the  lower  resonance  near  300  MHz.  The  gain  for  the 
horizontal  polarization  is  negligible  and  is  not  plotted.  The  gain  for  the  horizontal 


Figure  6.15:  Guided  power  of  a  microstrip  patch  antenna  over  a  uniform  and  over  a 
non-uniform  substrate. 
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Figure  6.16:  Radiation  pattern  of  a  microstrip  patch  antenna  over  a  uniform  and 
over  a  non-uniform  substrate  in  the  x-z  plane  for  the  vertical  polarization  at  the  first 
resonance. 

polarization  in  the  y-z  plane  is  plotted  in  Figure  6.17,  and  for  this  plane  the  vertical 
polarization  gain  is  negligible  and  not  plotted.  For  the  uniform  case,  the  gain  is  shown 
at  its  resonant  frequency  of  295  MHz,  and  for  the  non-uniform  case,  the  gain  is  shown 
at  its  resonant  frequency  of  300  MHz.  In  both  gain  plots,  the  non-uniform  case  has 
a  higher  gain  in  the  z  direction,  primarily  due  to  the  lower  loss  to  surface  waves. 

The  gain  for  the  vertical  polarization  in  the  x-z  plane  is  plotted  in  Figure  6.18 
for  both  configurations  at  the  second  resonance  which  is  near  410  MHz.  The  gain 
for  the  horizontal  polarization  is  negligible  and  is  not  plotted.  The  gain  in  the  y-z 
plane  is  plotted  in  Figure  6.19.  For  the  uniform  case,  the  gain  is  shown  at  its  resonant 
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Figure  6.17:  Radiation  pattern  of  a  microstrip  patch  antenna  over  a  uniform  and  over 
a  non-uniform  substrate  in  the  y-z  plane  for  the  horizontal  polarization  at  the  first 


resonance. 
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Figure  6.18:  Radiation  pattern  of  a  microstrip  patch  antenna  over  a  uniform  and  over 
a  non-uniform  substrate  in  the  x-z  plane  for  the  vertical  polarization  at  the  second 
resonance. 

frequency  of  409  MHz,  and  for  the  non-uniform  case,  the  gain  is  shown  at  its  resonant 
frequency  of  421  MHz.  In  both  gain  plots,  the  non-uniform  case  has  a  higher  gain  by 
about  a  factor  of  3,  due  to  the  significantly  lower  loss  to  surface  waves. 


GAIN 
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Figure  6.19:  Radiation  pattern  of  a  microstrip  patch  antenna  over  a  uniform  and  over 
a  non-uniform  substrate  in  the  y-z  plane  at  the  second  resonance. 
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6.4  Conclusions 

The  finite  difference-time  domain  technique  with  the  perfectly  matched  layer  ab¬ 
sorbing  boundary  conditions  successfully  analyzed  the  microstrip  patch  antenna  on 
a  non-uniform  substrate.  The  method  proved  to  be  flexible  and  accurate  in  han¬ 
dling  complex  dielectric  environments.  The  use  of  a  non-uniform  substrate  with  a 
microstrip  patch  antenna  has  been  shown  to  significantly  reduce  the  losses  due  to 
surface  wave  generation.  However,  the  addition  of  the  non-uniformity  can  change  the 
resonant  frequency  and  tends  to  reduce  the  bandwidth  of  the  antenna.  This  research 
only  begins  the  investigation  of  possible  uses  for  non-uniform  substrates  with  radiat¬ 
ing  elements.  Other  antennas  and  other  substrate  configurations  can  be  investigated 
with  only  minor  modifications  to  this  method. 


Chapter  7 
Conclusions 


The  finite  difference-time  domain  technique  has  been  used  to  successfully  analyze  a 
variety  of  electromagnetic  interference  and  radiation  problems.  Enhancements  to  the 
method  are  developed  and  implemented  which  increase  its  accuracy  and/or  efficiency. 
The  basic  FD-TD  algorithm  and  the  enhancements  which  were  developed  and  used 
are  described  in  detail.  The  technique  is  demonstrated  to  be  flexible  in  treating  this 
class  of  problems  which  include  potentially  complex  configurations  for  which  neither 
approximate  nor  analytical  methods  are  practical. 

All  the  electromagnetic  interference  problems  examined  concerned  potential  ra¬ 
diation  sources  within  a  computer  system.  The  design  and  implementation  of  com¬ 
ponents  within  computers  must  be  conducted  with  additional  consideration  for  elec¬ 
tromagnetic  emissions  because  regulatory  electromagnetic  emissions  limits  have  been 
extended  to  the  gigahertz  range  of  frequencies.  The  avoidance  of  potentially  efficient 
radiators  at  the  design  and  development  phases  is  of  interest  because  cost-prohibitive 
corrective  measures  , which  are  used  when  assembled  systems  encounter  problems,  can 


be  reduced. 
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Radiation  from  a  VLSI  heatsink  and  package  configuration  is  analyzed  using 
the  FD-TD  technique.  This  configuration  is  of  interest  for  the  following  reasons. 
Heatsinks  are  usually  metallic  and  in  close  proximity  to  a  chip  which  uses  large 
amounts  of  power.  The  heatsink  in  conjunction  with  the  ground  plane  will  have 
radiation  characteristics  similar  to  a  microstrip  patch  antenna.  The  electrical  dimen¬ 
sions  of  the  heatsink  are  comparable  to  the  wavelengths  of  concern.  In  light  of  the 
ever-increasing  power  levels  used  in  high-performance  integrated  circuits,  the  heatsink 
is  and  will  continue  to  be  an  indispensable  component  of  the  VLSI  chip  packaging 
configuration. 

The  heatsink  is  modeled  as  a  perfectly  conducting  rectangular  slab  positioned 
over  a  finite-size  dielectric  medium  representing  the  chip  package.  The  package  in 
turn  lies  on  a  dielectric  layer  on  a  ground  plane  of  infinite  extent,  which  models  the 
substrate  or  printed  wiring  board  with  at  least  one  reference  layer.  Simplified  excita¬ 
tion  sources  are  used  which  include  electric  and  magnetic  dipoles  oriented  vertically 
and  horizontally. 

In  applying  the  FD-TD  technique,  a  multi-zone  gridding  scheme  was  devel¬ 
oped  and  implemented.  The  multi-zone  gridding  scheme  implemented  allowed  a  high 
grid  resolution  in  the  vicinity  of  the  heatsink  without  sacrificing  overall  cell  number, 
thereby  allowing  increased  numerical  efficiency  for  low  frequencies.  This  discretization 
scheme  has  the  additional  advantages  of  improved  geometrical  and  material  modeling. 

The  analysis  of  the  heatsink,  a  critical  component  within  the  computer  packag¬ 
ing  environment,  illustrated  the  significance  of  resonance  due  to  appreciable  electrical 
dimensions  on  the  spurious  electromagnetic  radiation  from  the  VLSI  packaging  en- 


vironment.  Simulations  with  typical  heatsinks  dimensions  showed  the  occurrence  of 
resonance  in  the  low  gigahertz  range.  The  effects  of  the  presence  of  a  heatsink  on  the 
radiation  properties  of  dipole  models  have  been  explained  and  the  features  of  resonant 
behavior  have  been  presented.  The  effects  of  dielectric  loading  to  the  cavity  and  the 
addition  of  fins  to  the  heatsink  have  been  examined  and  quantified.  The  effective¬ 
ness  and  consequences  of  exploiting  practical  heatsink  implementation  options  such 
as  grounding  and  shielding  to  reduce  electromagnetic  emissions  have  been  discussed 
and  quantified  using  the  FD-TD  numerical  technique. 

The  radiation  properties  of  modules-on-backplane  configurations  have  been  an¬ 
alyzed  by  using  the  FD-TD  technique  and  by  taking  measurements.  The  modules- 
on-backplane  configuration  is  of  interest  because  it  is  commonly  used  in  computers. 
Also  the  reference  plane  dimensions  have  electrical  dimensions  which  are  appreciable 
at  frequencies  in  the  hundreds  of  megahertz  and  can  function  as  antennas. 

The  modules  and  backplanes,  which  are  printed  circuit  boards,  were  modeled  as 
perfectly  conducting  plates.  In  implementing  the  excitation  source,  a  voltage  source 
is  applied  between  the  module  and  the  backplane  at  the  location  of  the  connector. 
The  substrate  and  chips  populating  the  modules  were  not  considered  since  their  effect 
on  the  lowest  resonances  is  secondary  to  the  effect  on  the  resonances  from  the  ground 
planes  in  the  modules. 

An  extrapolation  scheme  based  on  Prony’s  method  was  used  in  conjunction  with 
a  relatively  short  FD-TD  time  response  to  efficiently  obtain  the  late  time  response. 
This  extrapolation  method  was  applied  to  the  time  response  of  the  current  and  voltage 
on  the  connector.  Prony’s  method  is  shown  to  be  useful  and  to  produce  accurate 
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results. 

Measurements  were  performed  using  an  HP  8510  Network  Analyzer.  The  net¬ 
work  analyzer  was  connected  to  the  modules-on-backplane  configuration  to  be  tested 
by  a  test  port  cable  and  semi-rigid  coaxial  cable.  In  each  of  the  measurements,  the 
center  conductor  of  the  coaxial  cable  was  soldered  to  the  module  under  test  such  that 
the  coaxial  cable  was  in  contact  with  at  least  one  of  the  reference  planes  in  the  mod¬ 
ule.  Fairly  good  agreement  was  obtained  between  the  measurements  and  the  FD-TD 
simulations. 

The  FD-TD  analysis  provided  radiation  impedances,  which  were  used  as  circuit 
model  parameters  to  estimate  overall  radiated  power.  The  analysis  of  the  modules-on- 
backplane  configuration  has  shown  that  radiation  enhancement  occurs  in  the  hundreds 
of  megahertz.  The  addition  of  modules  or  an  enclosure  have  been  shown  to  enhance 
the  radiation  from  the  structure.  The  potential  for  increased  radiation  due  to  the 
modules-on-backplane  configuration  has  been  demonstrated  with  measurements  and 
numerical  simulations. 

The  finite  difference-time  domain  technique  was  used  to  examine  the  coupling 
of  electromagnetic  energy  by  a  wire  penetrating  a  metallic  enclosure.  Metallic  en¬ 
closures  are  commonly  used  to  house  computer  systems  and  usually  have  cables  and 
wires,  which  transmit  data  or  power,  connected  to  them.  Since  the  metallic  enclosure 
functions  both  as  a  electromagnetic  shield  and  as  a  resonator  structure,  the  coupling 
of  energy  from  any  cables  exiting  from  the  box  is  especially  important.  The  currents 
on  the  interior  surface  of  the  enclosure  can  couple  to  wires  exiting  the  enclosure,  and 
the  wire  can  function  as  an  antenna  in  the  hundreds  of  MHz  frequency  range. 


Two  configurations  were  examined:  a  wire  extending  from  and  normal  to  a 
ground  plane  which  is  excited  by  a  coaxial  line,  and  a  grounded  wire  looping  up  from 
a  metal  screen  and  then  penetrating  the  metal  screen  and  extending  out  from  the 
screen  in  a  normal  direction  which  is  excited  by  a  dipole  on  the  side  of  the  screen  on 
which  the  wire  is  grounded.  Radiation  levels  associated  with  the  two  configurations 
were  quantified.  Significant  coupling  of  electromagnetic  energy  occured  when  the 
admittance  of  the  external  wire  was  matched  by  the  admittance  associated  with 
internal  side.  In  the  coaxial  line  configuration,  the  characteristic  admittance  of  the 
coaxial  line  matched  the  peaks  in  the  external  wire’s  radiation  admittance  well,  and 
resulted  in  radiated  power  only  1  to  2  dB  less  than  the  maximum  possible.  For  the 
infinite  metallic  screen  configurations,  fairly  significant  levels  of  radiation  occured 
at  frequencies  where  the  radiation  admittance  of  the  inner  wire  loop  matched  the 
radiation  admittance  of  the  external  wire  antenna.  For  the  two  cases  examined,  the 
maximum  coupling  was  only  3dB  and  6dB  less  than  the  power  radiated  by  the  dipole 
in  free  space. 

A  “smart”  cell  in  the  FD-TD  scheme,  which  can  model  the  radiation  admit¬ 
tance  of  a  wire  antenna,  has  been  developed  and  tested.  Analytical  solutions  for  the 
radiation  admittance  of  a  straight  wire  over  a  ground  plane  have  been  recited  and 
approximated  using  the  lumped  element  model,  which  consists  of  a  set  of  series  RLC 
circuits  in  parallel.  The  admittance  of  the  lumped  element  implemented  in  the  FD- 
TD  simulations  matched  the  circuit  solution  admittance  extremely  well.  Fairly  good 
correlation  between  lumped  element  FD-TD  simulations  and  full  FD-TD  simulations 
have  been  obtained  for  each  of  the  configurations  for  both  radiation  admittances  and 
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radiated  powers. 

One  clear  limitation  of  the  use  of  lumped  elements  in  modeling  radiation  admit¬ 
tances  is  that  the  results  were  only  as  accurate  as  the  lumped  element  model  of  the 
radiation  admittance.  The  lumped  element  model  considered  here  most  accurately 
modeled  the  peaks  in  the  radiation  conductance,  and  hence  produced  the  most  accu¬ 
rate  results  at  those  frequencies.  Correspondingly,  the  lumped  element  FD-TD  model 
had  increasing  errors  as  the  lumped  element  admittance  varied  from  the  theoretical 
radiation  admittance.  In  the  future,  if  methods  for  designing  passive  RLC  circuits 
to  more  accurately  approximate  arbitrary  admittances  as  a  function  of  frequency  are 
developed,  then  the  lumped  element  model  could  be  correspondingly  more  accurate. 

The  finite  difference-time  domain  technique  with  the  perfectly  matched  layer  ab¬ 
sorbing  boundary  conditions  successfully  analyzed  the  microstrip  patch  antenna  on 
a  non-uniform  substrate.  The  non-uniformity  considered  was  the  removal  of  a  rect¬ 
angular  ring  of  dielectric  from  the  substrate  which  surrounded  the  microstrip  patch 
antenna.  The  method  proved  to  be  very  flexible  in  handling  the  complicated  ge¬ 
ometries  considered  here.  The  perfectly  matched  layer  absorbing  boundary  condition 
displayed  excellent  absorbing  behavior  at  the  free  space  and  dielectric  interfaces. 

Radiation  impedance,  percentages  of  guided  power,  and  far-field  gain  patterns 
were  calculated,  presented  and  discussed.  The  use  of  a  non-uniform  substrate  with 
a  microstrip  patch  antenna  has  been  shown  to  significantly  reduce  the  losses  due  to 
surface  wave  generation.  However,  the  addition  of  the  non-uniformity  was  observed 
to  change  the  resonant  frequency  and  to  reduce  the  bandwidth  of  the  antenna. 


This  research  only  begins  the  investigation  of  possible  uses  for  non-uniform  sub¬ 
strates  with  radiating  elements.  Other  antennas  and  other  substrate  configurations 
can  be  investigated  with  only  minor  modifications  to  this  method.  In  particular,  the 
use  of  multiple  rings  instead  of  a  single  ring  of  free  space  could  provide  additional 
suppression  of  surface  waves.  An  alternative  to  non-uniform  substrates  would  be  peri¬ 
odic  substrates.  In  addition,  examining  arrays  of  antennas  and  their  mutual  coupling 
would  be  of  interest.  The  application  of  the  FD-TD  method  to  all  of  these  areas  is 
worthy  of  consideration. 

In  summary,  this  thesis  utilized  the  finite  difference-time  domain  technique  in 
analyzing  various  problems  and  discussed  enhancements  which  extended  the  method’s 
range  of  usefulness  by  increasing  its  accuracy  and/or  efficiency.  Important  issues  in 
electromagnetic  interference  in  computer  systems  were  examined  and  discussed.  The 
use  of  non-uniform  substrates  which  radiating  elements  was  examined  and  explored. 
Finally,  some  future  directions  for  research  in  these  areas  are  outlined. 
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